Transforming lat/lon coordinates into an undescribed DATUM from NOAA .grib2 files

Question

I seem to be having issues with converting coordinates.

I'm currently working in Python v3.7 osgeo (ogr, osr, gdal)

I want to find specific values (like temperature, wind, humidity) based on input coordinates in lat/lon (ex : 45N, -75E) that will be feed into a .grib2 provided by the NOAA weather forecast model called HRRR.

Any HRRR can be easily downloaded from : https://www.ftp.ncep.noaa.gov/data/nccf/com/hrrr/prod/hrrr.20200218/conus/ (I am specifically using the ones with the substring "hrrr.t00z.wrfsfcf")

Info from the HRRR file:

Driver: GRIB/GRIdded Binary (.grb, .grb2)

Files: hrrr.t00z.wrfsfcf00.grib2
       hrrr.t00z.wrfsfcf00.grib2.aux.xml

Size is 1799, 1059

Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",38.5],
    PARAMETER["standard_parallel_2",38.5],
    PARAMETER["latitude_of_origin",38.5],
    PARAMETER["central_meridian",262.5],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Metre",1]]

Origin = (-2699020.142521930392832,1588193.847443336388096)

Pixel Size = (3000.000000000000000,-3000.000000000000000)

Corner Coordinates:
Upper Left  (-2699020.143, 1588193.847) (134d 7'17.14"W, 47d50'44.64"N)
Lower Left  (-2699020.143,-1588806.153) (122d43'44.80"W, 21d 7'19.89"N)
Upper Right ( 2697979.857, 1588193.847) ( 60d53'28.48"W, 47d50'57.51"N)
Lower Right ( 2697979.857,-1588806.153) ( 72d16'48.48"W, 21d 7'28.62"N)
Center      (    -520.143,    -306.153) ( 97d30'21.52"W, 38d29'50.09"N)

I first begin by converting the lat/lon coordinates into the reference system of the file and then I try to find corresponding index.

from osgeo import ogr, osr, gdal

def get_reverse_point(coord_lat, coord_lon, in_epsg, output_srs_ref):
#From last comment in https://gis.stackexchange.com/questions/169224/getting-origin-as-long-lat-in-gdal-without-knowing-epsg
#This function converts lat/lon into the corresponding pixel coordinate of the file
#I also assume that the EPSG '4326' is the typical reference system used for lat/lon coordinates.


  in_epsg_ref = osr.SpatialReference()
  in_epsg_ref.ImportFromEPSG(in_epsg)

  coord_transform = osr.CoordinateTransformation(in_epsg_ref, out_srs_ref)

  # create a geometry from coordinates
  point = ogr.Geometry(ogr.wkbPoint)
  point.AddPoint(coord_lat, coord_lon) #Or : point.AddPoint(coord_lon, coord_lat)

  # transform point
  point.Transform(coord_transform)

  # return point in EPSG 4326
  return point.GetX(), point.GetY() #longitude and latitude?

Opened_HRRR_data = Open_HRRR_dataset(HRRR_file_full)
gt = Opened_HRRR_data.GetGeoTransform()

in_srs = osr.SpatialReference()
in_srs.ImportFromWkt(Opened_HRRR_data.GetProjection())

converted_points = get_reverse_point(45, -75, 4326, in_srs) 

print(converted_points)
>>>(1763296.6436729892, 940994.1197799345)

From the converted_points, which I assume gives me the corresponding pixels within the file, I try to compute its index with the following

def grib_data_get_index(grib_file, lat, lon, band_i=1):
    # From https://gis.stackexchange.com/questions/325924/getting-data-closest-to-specified-point-in-a-grib-file-python

    ox, pw, xskew, oy, yskew, ph = ds.GetGeoTransform()

    # calculate the indices (row and column)
    # Calculate the 4 nearest indices (i is for lats, j is for lons)
    i_bottom = math.floor((oy - lat) / ph)
    i_top = math.ceil((oy - lat) / ph)
    j_bottom = math.floor((lon - ox) / pw)
    j_top = math.ceil((lon - ox) / pw)

    # Save indices
    idx_array = np.array([i_bottom, i_top, j_bottom, j_top])

    # close the file
    del ds, band

    return idx_array

Which gives me the following array

HRRR_idx_array = grib_data_get_index(HRRR_file,converted_points[0],converted_points[1])
print(HRRR_idx_array)
>>> [ 58   59 1213 1214] #lat-lat-lon-lon equivalent

This doesn't seem so bad, and I assume so far I have placed the x/y, lat/lon inputs in their proper place (although I am ready to believe that I have not, but trying different combinations have only given me worst results), but when plotting the results on a map, it would seem as though the point is not really where it's suppose to be. The coordinate 45,-75 should be a bit more north-east to Lake Ontario, which is not the case looking at this plot. What could be wrong?

import matplotlib.pyplot as plt
band = Opened_HRRR_data.GetRasterBand(66) #band 66 is 2m Temperature
lc = band.ReadAsArray()
plt.imshow ( lc, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth )
plt.plot(HRRR_idx_array[2],HRRR_idx_array[0], 'ro') 

Answer 1

Ok so I was able to figure something out, but it's a completely different path than what I was trying to go with in this post.

First, I went with eccodes grib_ls. This allowed me to get the four nearest indices of a input lat/lon coordinate (my true original purpose). For example:

grib_ls -l 45.4,-74.32,1 $GRIB_FILE

This would give me four different values that would look like this:

• 1 - index=1552239 latitude=45.41 longitude=285.70 distance=1.59 (Km)
• 2 - index=1552238 latitude=45.41 longitude=285.66 distance=2.12 (Km)
• 3 - index=1550440 latitude=45.38 longitude=285.69 distance=2.22 (Km)
• 4 - index=1550439 latitude=45.39 longitude=285.65 distance=2.62 (Km)

Note that I don't know how to get eccodes to work in a Windows (which is what I had), maybe with cygwin, but I decided to install a windows subsystem for linux (WSL).

It would be too much of a hassle with eccodes to get the nearest values for each grib2 file, especially with what I wanted to have as a final result. However, if you're fluent with bash scripts then you don't need to follow the subsequent steps. In my case, Eccodes wasn't flexible enough, so I wanted to go back with my numpy array that has all the values of the variable in python. Now that I had the index, I needed to input that in my numpy array lc = band.ReadAsArray(). For that, I actually had to vertically flip the numpy array because the np array starts at the top left, while the eccodes index starts at the bottom left.

lc_flip = np.flip(lc,0)
lc_flip_flat = lc_flip.flatten()
print(lc_flip_flat[1552239])

And it works! Now to make it less manually dependent, it is possible with eccodes to save the output as a csv or a text file.

grib_ls -l 45.4,-74.32,1 $GRIB_FILE > tempo.txt

From this you just need to open the .csv or .txt in python and make it read the place where it mentions the indices and voilà!.