Chore/idsse 518/validate nbm temp

python/idsse_common/test/test_geo_image.py

@@ -21,7 +21,7 @@
 
 @fixture
 def proj():
-    proj_spec = '+proj=lcc +lat_0=25.0 +lon_0=-95.0 +lat_1=25.0 +r=6371200'
+    proj_spec = '+proj=lcc +lat_0=25.0 +lon_0=-95.0 +lat_1=25.0 +a=6371200'
     grid_spec = '+dx=2539.703 +dy=2539.703 +w=2345 +h=1597 +lat_ll=19.229 +lon_ll=-126.2766'
     return GridProj.from_proj_grid_spec(proj_spec, grid_spec)
 
@@ -382,19 +382,73 @@ def test_add_all_states(proj):
     # geo_image.show()
 
 
-# def test_temp_data(proj):
-#     colors = [(252, 180, 252), (187, 85, 192), (57, 21, 144), (0, 190, 242), (30, 186, 0),
-#               (254, 254, 0), (226, 46, 5), (149, 12, 6)]
-#     # colors = [(255, 0, 0), (0, 255, 0), (0, 0, 255)]
-#     current_path = os.path.dirname(os.path.realpath(__file__))
-#     filename = os.path.join(current_path, 'resources', 'nbm_temp-202211111100-202211121300.nc')
-#     attrs, data = read_netcdf(filename)
-#     print(attrs)
-#     if attrs['data_order'] == 'latitude,longitude':
-#         data = numpy.transpose(data)
-#     geo_image = GeoImage.from_data_grid(proj, data, ColorPalette.linear(colors),
-#                                         min_value=-40, max_value=80)
-#     geo_image.draw_state_boundary('All', color=(0, 0, 0))
-#     geo_image.show()
-
-#     assert False
+def test_temp_data(proj):
+    colors = [(252, 180, 252), (187, 85, 192), (57, 21, 144), (0, 190, 242), (30, 186, 0),
+              (254, 254, 0), (226, 46, 5), (149, 12, 6)]
+    # current_path = os.path.dirname(os.path.realpath(__file__))
+    # filename = os.path.join(current_path, 'resources', 'nbm_temp-202211111100-202211121300.nc')
+    filename = ('/Users/geary.j.layne/idssEngine/data/2024/01/12/'
+                'NBM.AWS.GRIB/TEMP/Fahrenheit/gridstore-74919423.nc')
+    # filename = ('/Users/geary.j.layne/idssEngine/data/2024/01/12/'
+    #             'NBM.AWS.GRIB/TEMP/Kelvin/gridstore-1518104402.nc')
+
+    attrs, data = read_netcdf(filename)
+    print(attrs)
+    if attrs['data_order'] == 'latitude,longitude':
+        data = numpy.transpose(data)
+    geo_image = GeoImage.from_data_grid(proj, data, ColorPalette.linear(colors),
+                                        min_value=-40, max_value=80)
+    geo_image.draw_state_boundary('All', color=(0, 0, 0))
+
+    # location of airport from https://airnav.com/airport/
+    locations = {
+        'KABQ': (-106.6082622, 35.0389316),
+        'KORD': (-87.9081497, 41.9769403),
+        'KMIA': (-80.2901158, 25.7953611)
+    }
+    # location of ASOS from https://www.faa.gov/air_traffic/weather/asos
+    locations = {
+        'KABQ': (-106.615, 35.042),
+        'KORD': (-87.932, 41.960),
+        'KMIA': (-80.317, 25.788)
+    }
+    # location of Sites from DESI (maybe via AWIPS)
+    locations = {
+        'KABQ': (-106.62, 35.05),
+        'KORD': (-87.93, 41.98),
+        'KMIA': (-80.28, 25.8)
+    }
+    for key, (lon, lat) in locations.items():
+        geo_image.outline_pixel(lon, lat, (0, 0, 0))
+        i, j = proj.map_geo_to_pixel(lon, lat, rounding='FLOOR')
+        geo_image.outline_pixel(i, j, (255, 255, 255), geo=False)
+        print(f'{key} ({lon}, {lat}) -> ({i}, {j}): {data[i, j]}')
+    geo_image.show()
+
+    # Convert Kelvin to Fahrenheit
+    # for i in range(864, 869):
+    #     for j in range(560, 565):
+    #         in_k = Q_(data[i, j], ureg.degK)
+    #         in_f = in_k.to('degF')
+    #         print(i, j, proj.map_pixel_to_geo(i, j),
+    #               data[i, j], in_k, in_f, in_f.to('degK'))
+
+    for (i, j) in [(865, 564), (1529, 865), (1867, 170)]:
+        print(i, j, ':', data[i, j])
+
+    # From DESI
+    # KABQ:
+    # lon, lat: -106.62, 35.05
+    # i, j: 865, 564
+    # 18 degF
+    # KORD:
+    # lon, lat: -87.93, 41.98
+    # i, j: 1529, 865
+    # 33 degF
+    # KMIA:
+    # lon, lat: -80.28, 25.8
+    # i, j: 1867, 170
+    # 76 degF
+
+    print(proj.map_geo_to_pixel(-105.2612, 39.9939))
+    assert False