[gdal-dev] handle GRIB files in polar stereographic projection

Fri Oct 24 05:37:59 PDT 2014

I was able to transform the grib files to NetCDF using the Unidata NetCDF
Tools. 
java -classpath netcdfAll-4.3.jar ucar.nc2.dataset.NetcdfDataset -in
ice_conc_nh.grb -out ice_conc_nh.nc

Amazingly GDAL can now work with this files. Here is the gdalinfo
(Coordinates are in Days)

$ gdalinfo NETCDF:ice_conc_nh.nc:Ice_cover_surface
Driver: netCDF/Network Common Data Format
Files: ice_conc_nh.nc
       ice_conc_nh.nc.aux.xml
Size is 760, 1120
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Stereographic"],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",-45.00000381469727],
    PARAMETER["scale_factor",0.9330127018922193],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["kilometre",1000,
        AUTHORITY["EPSG","9036"]]]
Origin = (-3698.443603515625000,6060.704589843750000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  Ice_cover_surface#coordinates=time 
  Ice_cover_surface#description=Ice cover
  Ice_cover_surface#Grib1_Center=88
  Ice_cover_surface#Grib1_Level_Desc=Ground or water surface
  Ice_cover_surface#Grib1_Level_Type=1
  Ice_cover_surface#Grib1_Parameter=91
  Ice_cover_surface#Grib1_Subcenter=0
  Ice_cover_surface#Grib1_TableVersion=1
  Ice_cover_surface#Grib_Variable_Id=VAR_88-0-1-91_L1
  Ice_cover_surface#grid_mapping=PolarStereographic_Projection
  Ice_cover_surface#long_name=Ice cover @ Ground or water surface
  Ice_cover_surface#missing_value=nan
  Ice_cover_surface#units=proportion
  NC_GLOBAL#Conventions=CF-1.6
  NC_GLOBAL#featureType=GRID
  NC_GLOBAL#file_format=GRIB-1
  NC_GLOBAL#GRIB_table_version=0,1
  NC_GLOBAL#history=Read using CDM IOSP GribCollection v2
  NC_GLOBAL#Originating_or_generating_Center=Oslo
  NC_GLOBAL#Originating_or_generating_Subcenter=0
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={1,6}
  NETCDF_DIM_time_VALUES=0
  PolarStereographic_Projection#earth_radius=6371229
  PolarStereographic_Projection#grid_mapping_name=stereographic
  PolarStereographic_Projection#latitude_of_projection_origin=90

PolarStereographic_Projection#longitude_of_projection_origin=-45.00000381469727

PolarStereographic_Projection#scale_factor_at_projection_origin=0.9330127018922193
  reftime#long_name=GRIB reference time
  reftime#standard_name=forecast_reference_time
  reftime#units=Hour since 2014-09-20T12:00:00Z
  time#long_name=GRIB forecast or observation time
  time#standard_name=time
  time#units=Hour since 2014-09-20T12:00:00Z
  x#standard_name=projection_x_coordinate
  x#units=km
  y#standard_name=projection_y_coordinate
  y#units=km
Corner Coordinates:
Upper Left  (   -3698.444,    6060.705) (166d23'34.54"E, 28d18'31.13"N)
Lower Left  (   -3698.444,   -5139.295) ( 80d44'25.12"W, 33d55'21.31"N)
Upper Right (    3901.556,    6060.705) (102d13'43.45"E, 27d32'42.72"N)
Lower Right (    3901.556,   -5139.295) (  7d47'44.26"W, 33d 1'11.09"N)
Center      (     101.556,     460.705) (122d34' 7.29"E, 85d27'18.90"N)
Band 1 Block=760x1 Type=Float32, ColorInterp=Undefined
  Min=-199.000 Max=100.000 
  Minimum=-199.000, Maximum=100.000, Mean=-51.973, StdDev=63.146
  NoData Value=nan
  Metadata:
    coordinates=time 
    description=Ice cover
    Grib1_Center=88
    Grib1_Level_Desc=Ground or water surface
    Grib1_Level_Type=1
    Grib1_Parameter=91
    Grib1_Subcenter=0
    Grib1_TableVersion=1
    Grib_Variable_Id=VAR_88-0-1-91_L1
    grid_mapping=PolarStereographic_Projection
    long_name=Ice cover @ Ground or water surface
    missing_value=nan
    NETCDF_DIM_time=0
    NETCDF_VARNAME=Ice_cover_surface
    STATISTICS_MAXIMUM=100
    STATISTICS_MEAN=-51.973378845782
    STATISTICS_MINIMUM=-199
    STATISTICS_STDDEV=63.146033919353
    units=proportion

This is the resulting Proj4-String
"+proj=stere +lat_0=90 +lon_0=-45.00000381469727 +k=0.9330127018922193
+x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=km +no_defs"

The current pitfall is that gdalwarp cannot transform the data

$ gdalwarp -t_srs EPSG:4326 NETCDF:ice_conc_nh.nc:Ice_cover_surface
trans/ice_conc_nh_4326.tiff
(image is attached)

Here is the gdalinfo of the WGS84 GeoTiff (SIZE gets very small &
Coordinates are in Days like the original data)

Driver: GTiff/GeoTIFF
Files: trans/ice_conc_nh_4326.tiff
Size is 2761, 217
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-176.400770410834753,55.428082489041401)
Pixel Size = (0.128740639134764,-0.128740639134764)
... ... ...Metadata... ... ...
Corner Coordinates:
Upper Left  (-176.4007704,  55.4280825) (176d24' 2.77"W, 55d25'41.10"N)
Lower Left  (-176.4007704,  27.4913638) (176d24' 2.77"W, 27d29'28.91"N)
Upper Right ( 179.0521342,  55.4280825) (179d 3' 7.68"E, 55d25'41.10"N)
Lower Right ( 179.0521342,  27.4913638) (179d 3' 7.68"E, 27d29'28.91"N)
Center      (   1.3256819,  41.4597231) (  1d19'32.45"E, 41d27'35.00"N)
... ... ...Other Stuff... ... ...

Do you see a way to fix the data and transform it correctly?

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/handle-GRIB-files-in-polar-stereographic-projection-tp5166227p5169393.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.