diff --git a/python/idsse_common/test/test_geo_image.py b/python/idsse_common/test/test_geo_image.py index 5fa8da4e..51d5b6ce 100644 --- a/python/idsse_common/test/test_geo_image.py +++ b/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