5

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') 

2m temperature map with the red point pointing at the supposed 45N,-75E coordinate (not accurate)]

2

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à!.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.