Gdal grib import Lambert problem, gets origin and corners wrong

[alexice]

Dear developers,

When importing a grib file with gdal to grass that contains a Lambert Projection, currently the Latitude of Origin is set to zero despite the fact that the grib file actually contains this info.

Therefore I suggest to change the code around line # 487 in gribdataset.cpp to:

      case GS3_LAMBERT:
        oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
                    meta->gds.meshLat, meta->gds.orientLon,
                    0.0, 0.0); // set projection
        break;

now gds.meshLat is used to set the Latitude of Origin within the projection.

Another problem occurs when the grib import routine is reconstructing the corners of the domain. In line # 559 of gribdataset.cpp the calculation is correct to go from the lower left to the upper left corner:

           rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel

After calculating the upper left corner, the grib input sends everything off to gdal and gdal calculates the lower left corner using the just calculated upper left corner and nRasterYSize. Now this causes a problem, because we just used gds.Ny -1 steps to go up to the upper left corner, but used gds.Ny steps to go back to the lower left corner. nRasterYSize is defined in line # 465 as

    nRasterYSize = meta->gds.Ny;

Same problem exists in the x direction.

As an example one can use the grib file AWIP32-fixed.grb from the NARR dataset, which contains a Lambert conformal projection. The file can be found under: ​http://nomads.ncdc.noaa.gov/data/narr/AWIP32-fixed.grb. Using this file, gds.Ny is 277 and gds.Nx is 349. So currently 276 steps are used to go to the upper left corner of the domain (which is correct, given 277 points) and 277 steps are used to reconstruct the lower left corner.

Running gdalinfo on the AWIP32-fixed file shows what happens.

$ gdalinfo AWIP32-fixed.grb
Driver: GRIB/GRIdded Binary (.grb)
Files: AWIP32-fixed.grb
Size is 349, 277
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371200,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",50],
    PARAMETER["standard_parallel_2",50],
    PARAMETER["latitude_of_origin",50],
    PARAMETER["central_meridian",-107],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-5632642.225474951788783,4347242.348625713959336)
Pixel Size = (32463.000000000000000,-32463.000000000000000)
Corner Coordinates:
Upper Left  (-5632642.225, 4347242.349) (148d38'24.25"E, 46d38'4.34"N)
Lower Left  (-5632642.225,-4645008.651) (145d23'44.80"W,  0d48'55.65"N)
Upper Right ( 5696944.775, 4347242.349) (  2d29'41.57"W, 46d 3'58.26"N)
Lower Right ( 5696944.775,-4645008.651) ( 68d14'23.82"W,  0d36'32.18"N)
Center      (   32151.275, -148883.151) (106d33'44.57"W, 48d39'37.36"N)

The lower left corner is put to 145d23'44.80"W,0d48'55.65"N whereas the definition of the lower left corner was at -145.5W 1N as shown here in the output of degrib:

GDS | Number of Points | 96673
GDS | Projection Type | 3 (Lambert Conformal)
GDS | Shape of Earth | sphere
GDS | Radius | 6371.200000 (km)
GDS | Nx (Number of points on parallel) | 349
GDS | Ny (Number of points on meridian) | 277
GDS | Lat1 | 1.000000
GDS | Lon1 | -145.500000
GDS | u/v vectors relative to | grid
GDS | Dx | 32463.000000 (m)
GDS | Dy | 32463.000000 (m)
GDS | Input GRIB2 grid, scan mode | 64 (0100)
GDS | Output grid, scan mode | 64 (0100)
GDS | (.flt file grid), scan mode | 0 (0000)
GDS | Output grid, scan i/x direction | positive
GDS | Output grid, scan j/y direction | positive
GDS | (.flt file grid), scan j/y direction | negative
GDS | Output grid, consecutive points in | i/x direction
GDS | Output grid, adjacent rows scan in | same direction
GDS | MeshLat | 50.000000
GDS | OrientLon | -107.000000
GDS | Which pole is on the plane | North
GDS | bi-polar projection | No
GDS | Tangent Lat1 | 50.000000
GDS | Tangent Lat2 | 50.000000
GDS | Southern Lat | 0.000000
GDS | Southern Lon | 0.000000

Can this problem be fixed easily? Let me know if I can help.

Cheers, Alexice

[Jukka Rahkonen]

The grib file ​http://nomads.ncdc.noaa.gov/data/narr/AWIP32-fixed.grb is still available and gdalinfo (GDAL 2.0-dev) shows still the same report with wrong corner coordinates.

[Kyle Shannon]

I don't think this is wrong. The GRIB origin is considered the center of the lower left pixel. It's used as the base coordinate for projected origins. The GDAL lower left is 1/2 pixel in the -x and -y directions: which give the coordinates you are seeing:

reported corner:

kyle@kyle-workstation:~/Desktop/gdal/2607$ gdaltransform -s_srs '+proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs ' -t_srs EPSG:4326
-5648873.725 -4628777.151
-145.540043353481 0.855548456209198 0

reported corner + one half cell in x and y:

kyle@kyle-workstation:~/Desktop/gdal/2607$ gdaltransform -s_srs '+proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs ' -t_srs EPSG:4326
-5565943.517711081, -4649977.057317405
-145 0.999999999999988 0

Note that the coordinates I used differ from yours, but are from gdalinfo in the trunk.

It's a difference in data model, and IMHO not a bug. Inclined to close as wontfix.

Please reopen if needed.

[Jukka Rahkonen]

Embarrassing, did not pay attention to "W" in 145d23'44.80"W and thought it had wrong sign...