1

I'm trying to work with NOAA's NDFD weather data, which are released as rasters in GRIB format. Specifically, I'm currently trying to get the pixels that intersect an arbitrary geometry.

Information about the grid's projection is here, and an example GRIB file is here (42MB).

When I load the file with rasterio, I see that the bounds are essentially 0-1377 in the y direction and 0-2145 in the x direction:

import rasterio
r = rasterio.open('YEUZ98_KWBN_201701010519')
r.bounds
>>> BoundingBox(left=-0.5, bottom=1376.5, right=2144.5, top=-0.5)

When I take that same file and run gdalsrsinfo on it, I get a proj4 string:

> gdalsrsinfo YEUZ98_KWBN_201701010519 -o proj4

+proj=lcc +lat_0=25 +lon_0=265 +lat_1=25 +lat_2=25 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs

However if I take an arbitrary point in the US, and try to project it to that projection, I get x and y values way outside the 0-1377, 0-2145 range I expect:

import geopandas as gpd
import pyproj
from shapely.geometry import Point

proj_str = '+proj=lcc +lat_0=25 +lon_0=265 +lat_1=25 +lat_2=25 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs'
crs = pyproj.CRS.from_proj4(proj_str)

point = Point(-122.4194, 37.7749)
gdf = gpd.GeoDataFrame(geometry=[point], crs={'init': 'epsg:4326'})
gdf = gdf.to_crs(crs)
list(gdf.iloc[0].geometry.coords)

# [(-2456716.706510244, 1682088.477346224)]

I expected this code to reproject the WGS84 point to a point on the NDFD grid, but it doesn't. Why not?

Edit: output of gdalinfo:

> gdalinfo YEUZ98_KWBN_201701010519
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: YEUZ98_KWBN_201701010519
Size is 2145, 1377
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Coordinate System imported from GRIB file",
        DATUM["unnamed",
            ELLIPSOID["Sphere",6371200,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",265,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["Metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["Metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-2764474.350731994025409,3232111.710792394354939)
Pixel Size = (2539.702999999999975,-2539.702999999999975)
Corner Coordinates:
Upper Left  (-2764474.351, 3232111.711) (130d 7'22.60"W, 49d56'49.26"N)
Lower Left  (-2764474.351, -265059.320) (121d33'48.72"W, 20d10'43.04"N)
Upper Right ( 2683188.584, 3232111.711) ( 60d51'57.95"W, 50d 6'46.52"N)
Lower Right ( 2683188.584, -265059.320) ( 69d11'54.73"W, 20d19' 6.40"N)
Center      (  -40642.883, 1483526.195) ( 95d27' 8.65"W, 38d13' 5.87"N)
Improve this question
9
  • "the bounds are essentially 0-1377 in the y direction and 0-2145 in the x direction". The example is 2145 columns x 1377 rows, so rasterio isn't picking up the georeferencing and is just giving you the pixel coordinates, not the georeferenced coordinates. What is the output of gdalinfo examplegrib? gdalinfo tells me it is properly georeferenced.
    – user2856
    Commented Dec 21, 2019 at 5:24
  • Yeah, I was also confused by that, because I thought rasterio just wraps GDAL, but gdalinfo says it's georeferenced and rasterio doesn't? I've been searching for a way to manually pass a crs to rasterio when reading a file but haven't found it yet.
    – Kyle Barron
    Commented Dec 21, 2019 at 5:32
  • But regardless of the rasterio issue, (which might be this issue), there must exist some transformation that takes latitude and longitude WGS84 points and converts them to the grid defined by these parameters, right?
    – Kyle Barron
    Commented Dec 21, 2019 at 5:41
  • Use with rasterio.Env(), rasterio.read('YEUZ98_KWBN_201701010519') as r:
    – user2856
    Commented Dec 21, 2019 at 5:45
  • And you could do your point coord reprojection in rasterio to rasterio.readthedocs.io/en/stable/api/…
    – user2856
    Commented Dec 21, 2019 at 6:16

1 Answer 1

Reset to default
3

I expected this code to reproject the WGS84 point to a point on the NDFD grid, but it doesn't. Why not?

Your code does reproject the WGS84 point to a point on the NDFD grid.

Your problem is that rasterio is not picking up your GRIB file georeferencing and just returning the bounds as pixel coordinates from 0-ncols, 0-nrows (in fact your code should not work at all as there is no rasterio.read function, but I'll assume you meant rasterio.open).

with rasterio.open('YEUZ98_KWBN_201701010519') as r:

    print(r.bounds)
    # BoundingBox(left=-0.5, bottom=1376.5, right=2144.5, top=-0.5)           

    print(r.shape)
    # (1377, 2145)

    print(r.crs)
    # None

You need to ensure your GDAL_DATA environment is set correctly (either as an environment variable that is picked up implicitly or explicitly via with rasterio.Env(GDAL_DATA='/path/to/gdal/data')) and then wrap your call to rasterio.open with the rasterio.Env() context manager.

Then you can reproject your points to the NDFD grid (either using your shapely and geopandas method or various other ways including rasterio.warp.transform) and intersect your GRIB with your reprojected points. For example:

import rasterio
from rasterio.warp import transform

with rasterio.Env(), rasterio.open('YEUZ98_KWBN_201701010519') as r:
## OR with rasterio.Env(GDAL_DATA='/path/to/gdal/data'), rasterio.open('YEUZ98_KWBN_201701010519') as r:

    lon_lats = [[-122.4194, 37.7749], [-122.5, 37.75]]
    xys = list(zip(*transform({'init': 'epsg:4326'}, r.crs, *zip(*lon_lats))))
    data = list(r.sample(xys))

    print(r.bounds)
    # >>> BoundingBox(left=-2764474.3507319884, bottom=-265059.3202076056, right=2683188.584268011, top=3232111.7107923944)

    print(xys)
    # >>> [(-2456716.706510244, 1682088.477346224), (-2464411.5901549244, 1680767.699278633)]

    print(data)
    # >>> [array([ 8.85      ,  8.35      ,  8.35      ,  7.7499939 ,  7.7499939 ,
    #          7.7499939 ,  7.7499939 ,  7.7499939 ,  7.2499939 ,  7.2499939 ,
    #          7.2499939 ,  8.35      ,  8.85      ,  9.9500061 , 10.55001221,
    #         11.15001831, 11.65001831, 11.65001831, 11.65001831, 10.55001221,
    #          9.4500061 ,  8.85      ,  8.85      ,  8.35      ,  7.7499939 ,
    #          7.7499939 ,  7.2499939 ,  7.2499939 ,  6.65001831,  6.65001831,
    #          6.65001831,  6.65001831,  6.15001831,  6.15001831,  6.15001831,
    #          7.2499939 ,  7.7499939 ,  9.4500061 ,  9.9500061 ,  8.85      ,
    #          8.35      ,  7.7499939 ,  7.2499939 ,  7.2499939 ,  8.35      ,
    #          9.9500061 , 10.55001221]),
    #  array([ 7.7499939 ,  7.2499939 ,  7.7499939 ,  7.2499939 ,  7.2499939 ,
    #          6.65001831,  6.65001831,  6.65001831,  6.15001831,  6.15001831,
    #          6.15001831,  7.7499939 ,  8.85      ,  9.4500061 , 10.55001221,
    #         11.15001831, 11.15001831, 11.15001831, 10.55001221,  9.4500061 ,
    #          8.35      ,  8.35      ,  7.7499939 ,  7.2499939 ,  7.2499939 ,
    #          6.65001831,  7.2499939 ,  6.65001831,  6.65001831,  6.15001831,
    #          6.15001831,  6.15001831,  5.55001221,  5.55001221,  5.55001221,
    #          6.65001831,  7.7499939 ,  9.9500061 ,  9.4500061 ,  7.7499939 ,
    #          6.65001831,  6.15001831,  6.15001831,  5.55001221,  8.35      ,
    #         10.55001221,  9.9500061 ])]

If you still can't get rasterio to read the georeferencing, wrap the GRIB in a VRT:

gdal_translate -of VRT YEUZ98_KWBN_201701010519 YEUZ98_KWBN_201701010519.vrt

And rasterio will pick up the georeferencing without the Env context manager:

with rasterio.open('YEUZ98_KWBN_201701010519.vrt') as r:

    print(r.crs.to_proj4())
    # >>> +proj=lcc +lat_1=25 +lat_2=25 +lat_0=25 +lon_0=265 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs=True

    print(r.bounds)
    # >>> BoundingBox(left=-2764474.3507319884, bottom=-265059.3202076056, right=2683188.584268011, top=3232111.7107923944)
Improve this answer
2
  • Thanks! After a few more hours working with those last comments, I was able to figure it out. One thing that was confusing me was that even after using rasterio.Env, the projected units were in meters and not in cells.
    – Kyle Barron
    Commented Dec 23, 2019 at 0:31
  • @KyleBarron, not sure what you mean. The units are supposed to be metres, not cells. There are functions in rasterio to convert pixel/cell coordinates to map/georeferenced coordinates and vice-versa if you need that? But the sample function handles that for you when extracting pixel data values.
    – user2856
    Commented Dec 23, 2019 at 2:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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