Opened 10 years ago

Closed 8 years ago

#4807 closed defect (fixed)

gdalinfo reports wrong latitude of origin for common GRIB2 lambert conformal rasters

Reported by: kempisty Owned by: warmerdam
Priority: normal Milestone: 1.10.0
Component: default Version: 1.9.0
Severity: normal Keywords: gdalinfo grib grib2 raster proj4 lambert
Cc: Kyle Shannon

Description

The National Weather Service produces forecasts for the contiguous United States (CONUS) on a Lambert conformal raster, and delivers them in GRIB2 format ( http://www.nws.noaa.gov/mdl/degrib/dataloc.php#ndfd )

"gdalinfo -proj4" reports "+lat_0=0" for these files, which is incorrect. It should read "+lat_0=25".

PROJCS string also reports "latitude_of_origin=0", which should be 25.

Here is the GRIB2 spec: http://graphical.weather.gov/docs/grib_design.html

The correct lat_0 / latitude_of_origin for a GRIB2 Lambert raster is found in "Section 3 - Grid Definition Section", octet 48-51.

Change History (6)

comment:1 by Even Rouault, 10 years ago

Did you notice an *actual* issue with the georeferencing of those datasets ? I've downloaded http://weather.noaa.gov/pub/SL.us008001/ST.opnl/DF.gr2/DC.ndfd/AR.conus/VP.001-003/ds.apt.bin and as lat_1 = lat_2 (= 25), I believe that the value of lat_0 is ignored. I've tested that by assigning lat_0 with the meshLat field (whose value is 25), and I get the same georeferencing when looking at the coordinates in latitude,longitude on gdalinfo output. So this might be more correct, but I'm a bit hesitant to apply not being neither a GRIB specialist nor a LCC one...

Index: frmts/grib/gribdataset.cpp
===================================================================
--- frmts/grib/gribdataset.cpp	(révision 24921)
+++ frmts/grib/gribdataset.cpp	(copie de travail)
@@ -726,7 +726,7 @@
         break;
       case GS3_LAMBERT:
         oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
-                    0.0, meta->gds.orientLon,
+                    meta->gds.meshLat, meta->gds.orientLon,
                     0.0, 0.0); // set projection
         break;

comment:2 by kempisty, 10 years ago

Yes, unfortunately it does cause problems. The problems appear when warping the rasters to other map projections. I actually noticed that there was a problem when I used gdalwarp to convert these rasters to Spherical Mercator, for use with Google Maps. The resulting mercator raster should have had a perfectly horizontal boundary across the northern US, where the US/Canada border runs along 49N. Instead of horizontal, that boundary was both curved and sloped. That's when I said, "Something is wrong here," and I traced it to the lat_0 value.

Once I ignored gdalinfo's lat_0=0 and changed it to lat_0=25, I had the mercator I was looking for.

Our office should probably do more to help maintain the grib2 driver, frankly. I will bring it up with the powers that be, but you know how the gov budget is... also curved and sloped, and in the wrong direction. :-(

comment:3 by Even Rouault, 10 years ago

You're right. I've just tested warping and indeed, the result is OK when lat_0=25, but not when lat_0=0. What is weird is that the extent boundary reported by gdalinfo are the same in both cases, but oh well...

So I've applied the above patch in trunk in r24928. As I didn't test with other GRIB datasets (possibly from other sources than the National Weather Service), I didn't apply in 1.9 branch. Perhaps if you have the opportunity to test, you could report if that works OK with other datasets too

comment:4 by Kyle Shannon, 8 years ago

Cc: Kyle Shannon added

comment:5 by Jukka Rahkonen, 8 years ago

No bad feedback about r24928 and I suppose that it is now included in recent GDAL versions so probably ticket is mature enough for closing. Can anybody (Even?) confirm?

comment:6 by Even Rouault, 8 years ago

Milestone: 1.10.0
Resolution: fixed
Status: newclosed
Note: See TracTickets for help on using tickets.