Spatial reference system transformation differs from wgrib2

[jswhit]

From Cheewai....@gmail.com on May 07, 2010 02:18:12

I have a GRIB2 data file ( ftp://ftp.ssec.wisc.edu/pub/eosdb/dbcras/2010_05_06_126/06May10.12z.dbCRAS G.SouthAfrica.grib2) which contains grid data (140x210) in polar stereographic projection.

It appears that pygrib latlons() method automatically reprojects the grid point coordinates from stereographic to geographic (latitude/longitude).

The 'wgrib2' utility ( http://www.cpc.noaa.gov/products/wesley/wgrib2/ ) is also capable of transforming the grid point coordinates.

But, some return values of latlons() differ from those in the CSV dump of 'wgrib2'.

For ease of comparison, I have put the coordinate values side-by-side in a single text file (dbcras-lonlats-compare.txt) in the following format:

"lon-lat from wgrib2" "lon-lat from pygrib"

To illustrate, below is a snippet the file: -4.50056,-45.7956 -4.514481,-45.789450 -4.02109,-46.0063 -4.035438,-46.000016 -3.53655,-46.2144 -3.551329,-46.208080 -3.04692,-46.4199 -3.062145,-46.413584 -2.55221,-46.6228 -2.567878,-46.616468 -2.0524,-46.8231 -2.068522,-46.816676 -1.54749,-47.0206 -1.564075,-47.014148 -1.03749,-47.2153 -1.054539,-47.208826 -0.522396,-47.4071 -0.539918,-47.400650 -0.00222074,-47.596 -0.020220,-47.589561

As you can see, some of them vary quite significantly.

I am a newbie in Python and GIS field. Do you know what could possibly account for the variation?

I am using 32-bit Ubuntu 10.04 with wgrib2 v0.1.8.2 and pygrib SVN revision 192 .

Attachment: dbcras-lonlats-compare.txt

Original issue: http://code.google.com/p/pygrib/issues/detail?id=6

[jswhit]

From whitaker.jeffrey@gmail.com on May 07, 2010 05:00:30

The latlons method always returns the lats and lons (as the name suggests). If you want the projection coordinates, you can do

import pyproj P = pyproj.Proj(grb.projparams) # grb is the grib message object x, y = P(lons, lats)

As for why the result differs in wgrib2, without digging through the wgrib2 code I have no idea. I will try to investigate when I get some time.

-Jeff

[jswhit]

From whitaker.jeffrey@gmail.com on May 07, 2010 08:05:16

wgrib2 is using it's own custom code for map projections. pygrib uses the industry standard, battle tested proj4 library. In other words, I trust pygrib's answer more than wgrib2.

However, if you have evidence that wgrib2's answer is correct, and there is a bug in pygrib, I will investigate further.

-Jeff

[jswhit]

From whitaker.jeffrey@gmail.com on May 27, 2010 08:35:26

Here's the response from the wgrib2 developer Wes Ebisuzaki:

Jeff,

I just looked at the difference between the lat-lon from wgrib2 and pygrib in https://code.google.com/p/pygrib/issues/detail?id=6 Since the differences are small, I suspected that the codes were using a different radius of the earth

and here is what I found.

(1) wgrib2 used the default earth radius in the case where the user specified a radius fixed - will be in the next release

(2) the grib file specified that the radius of the earth was 637.1189 km. (scaling factor was off)

(3) When I changed the radius of the earth to 6371.189 km, wgrib2 got the same results

(4) you should check the pygrib code to see why it got the correct radius even though the radius was encoded incorrectly.

             Wesley

[jswhit]

From whitaker.jeffrey@gmail.com on May 27, 2010 12:34:35

Wesley is right, and there is indeed a bug in the grib file. The scale factor for the specified earth radius is set to 1/10, but the unscaled radius is provided in meters. pygrib was getting the right answer for the wrong reason by ignoring the scale factor. This is really a bug in pygrib though, and if I fix it you will get nonsensical values from the latlons method (since the earth radius will be set to 637 km).

[jswhit]

From Cheewai....@gmail.com on May 27, 2010 20:56:40

Thank you all for your effort.

Jeff, do fix pygrib. It's better to bite the bullet and do the "right" thing. ;-)

I will forward your findings to the dbCRAS folks at SSEC who produced the grib file.

[jswhit]

From whitaker.jeffrey@gmail.com on May 28, 2010 04:55:06

I have fixed the bug, but added a workaround that forces the earth radius to be O(10**6) meters. With this,the latlons method still produces the right answer with the dbcras files.

[jswhit]

From fepeacoc...@googlemail.com on May 06, 2012 04:36:11

Hello Jeff

Continuing on the same theme: I compared wgrib2's v.1.9.5.1 grid coordinates for EUMETSAT's MPE grib2 file with that using pygrib's v1.9.4 and found similar inconsistencies. Is it true to say that pygrib's coordinates are correct since wgrib2 is still in ALPHA stage and pygrib had the functionality coded since v1.5 and based on the previous comments the same problem may exist in wgrib2. Is it correct to assume that the earth radius may play a part in this issue as well and that wgrib2 has not improved after the 1.9.5.1 update?

[jswhit]

From whitaker.jeffrey@gmail.com on May 06, 2012 04:54:55

I don't know if Wesley has updated wgrib2. I think the above discussion speaks for itself - the fundamental problem is that the grib file is miscoded. pygrib2 corrects for this, wgrib2 does not.

[jswhit]

From whitaker.jeffrey@gmail.com on December 29, 2014 11:57:22

Status: Done