Opened 15 years ago
Closed 6 years ago
#2584 closed defect (fixed)
NetCDF Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover flipped east/west?
Reported by: | Markus Neteler | Owned by: | warmerdam |
---|---|---|---|
Priority: | normal | Milestone: | |
Component: | GDAL_Raster | Version: | svn-trunk |
Severity: | normal | Keywords: | netcdf |
Cc: | dnadeau |
Description
Importing maps from leads to the funny effect that east/west are flipped.
Note that some funky characters are reported by gdalinfo below.
Procedure:
wget ftp://sidads.colorado.edu/pub/DATASETS/snow/nsidc0321v01/north/NL.2004.nsidc0321v01.nc gdalinfo NL.2004.nsidc0321v01.nc Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute Driver: netCDF/Network Common Data Format Files: NL.2004.nsidc0321v01.nc Size is 512, 512 Coordinate System is `' Metadata: NC_GLOBAL#creation_date=Wed Feb 20 13:19:33 MST 2008 NC_GLOBAL#grid_id=Nl NC_GLOBAL#references=Data set documentation: http://nsidc.org/data/nsidc-0321.html NC_GLOBAL#comment=Not set at this time NC_GLOBAL#source=Satellite observations from DMSP SSM/I and Terra MODIS NC_GLOBAL#institution=National Snow & Ice Data Center, Boulder, CO NC_GLOBAL#title=Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover NC_GLOBAL#history=Wed Feb 20 13:20:20 2008: ncrcat -O NL.2004001.nsidc0321v01.nc NL.2004009.nsidc0321v01.nc NL.2004017.nsidc0321v01.nc NL.2004025.nsidc0321v01.nc NL.2004033.nsidc0321v01.nc NL.2004041.nsidc0321v01.nc NL.2004049.nsidc0321v01.nc NL.2004057.nsidc0321v01.nc NL.2004065.nsidc0321v01.nc NL.2004073.nsidc0321v01.nc NL.2004081.nsidc0321v01.nc NL.2004089.nsidc0321v01.nc NL.2004097.nsidc0321v01.nc NL.2004105.nsidc0321v01.nc NL.2004113.nsidc0321v01.nc NL.2004121.nsidc0321v01.nc NL.2004129.nsidc0321v01.nc NL.2004137.nsidc0321v01.nc NL.2004145.nsidc0321v01.nc NL.2004153.nsidc0321v01.nc NL.2004161.nsidc0321v01.nc NL.2004169.nsidc0321v01.nc NL.2004177.nsidc0321v01.nc NL.2004185.nsidc0321v01.nc NL.2004193.nsidc0321v01.nc NL.2004201.nsidc0321v01.nc NL.2004209.nsidc0321v01.nc NL.2004217.nsidc0321v01.nc NL.2004225.nsidc0321v01.nc NL.2004233.nsidc0321v01.nc NL.2004241.nsidc0321v01.nc NL.2004249.nsidc0321v01.nc NL.2004257.nsidc0321v01.nc NL.2004265.nsidc0321v01.nc NL.2004273.nsidc0321v01.nc NL.2004281.nsidc0321v01.nc NL.2004289.nsidc0321v01.nc NL.2004297.nsidc0321v01.nc NL.2004305.nsidc0321v01.nc NL.2004313.nsidc0321v01.nc NL.2004321.nsidc0321v01.nc NL.2004329.nsidc0321v01.nc NL.2004337.nsidc0321v01.nc NL.2004345.nsidc0321v01.nc NL.2004353.nsidc0321v01.nc NL.2004361.nsidc0321v01.nc NL.2004.nsidc0321v01.nc NC_GLOBAL#nco_openmp_thread_number=1 Subdatasets: SUBDATASET_1_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":bpInfo SUBDATASET_1_DESC=[46x81] bpInfo (XÄ���ر����enmp_thread_num���) SUBDATASET_2_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":scaInfo SUBDATASET_2_DESC=[46x81] scaInfo (XÄ���ر����enmp_thread_num���) SUBDATASET_3_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":latitude SUBDATASET_3_DESC=[721x721] latitude (64-bit floating-point) SUBDATASET_4_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":longitude SUBDATASET_4_DESC=[721x721] longitude (64-bit floating-point) SUBDATASET_5_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":SWE SUBDATASET_5_DESC=[46x721x721] lwe_thickness_of_surface_snow_amount (8-bit integer) SUBDATASET_6_NAME=NETCDF:"NL.2004.nsidc0321v01.nc":SCA SUBDATASET_6_DESC=[46x721x721] surface_snow_area_fraction (8-bit integer) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0) Warning 1: Unable to save auxilary information in NL.2004.nsidc0321v01.nc.aux.xml. # 8-day Snow-Covered Area (SCA) gdal_translate NETCDF:"NL.2004.nsidc0321v01.nc":SCA /tmp/sca.tif qgis /tmp/sca.tif -> see attachment
The projection is described here: http://nsidc.org/data/ease/ease_grid.htm (and looks better than what QGIS sees).
A bug?
Markus
Attachments (3)
Change History (17)
by , 15 years ago
Attachment: | gdal_nc_ease_grid_bug.png added |
---|
comment:1 by , 15 years ago
Cc: | added |
---|
by , 15 years ago
Attachment: | netcdf_flipped_by_gdal.png added |
---|
by , 15 years ago
Attachment: | netcdf_flipped_by_gdal.2.png added |
---|
ETOPO1 netCDF Grid is vertical-flipped by GDAL
comment:2 by , 15 years ago
I have a similar problem, with ETOPO1 netCDF grid.
Procedure:
wget http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/ETOPO1_Ice_g.grd.gz gunzip ETOPO1_Ice_g.grd.gz gdalinfo ETOPO1_Ice_g.grd Driver: netCDF/Network Common Data Format Files: ETOPO1_Ice_g.grd Size is 21601, 10801 Coordinate System is `' Origin = (-180.008333333333326,90.008333333333340) Pixel Size = (0.016666666666667,-0.016666666666667) Metadata: NC_GLOBAL#Conventions=COARDS/CF-1.0 NC_GLOBAL#title=ETOPO1_Ice_g.int.grd NC_GLOBAL#history=grdreformat ETOPO1_Ice_g.int.grd=ni ETOPO1_Ice_g.grd=ns NC_GLOBAL#GMT_version=4.3.1b NC_GLOBAL#node_offset=0 z#long_name=z z#_FillValue=-32768 z#actual_range=-10898, 8271 lon#long_name=longitude lon#units=degrees_east lon#actual_range=-180, 180 lat#long_name=latitude lat#units=degrees_north lat#actual_range=-90, 90 Corner Coordinates: Upper Left (-180.0083333, 90.0083333) Lower Left (-180.0083333, -90.0083333) Upper Right ( 180.0083333, 90.0083333) Lower Right ( 180.0083333, -90.0083333) Center ( 0.0000000, 0.0000000) Band 1 Block=21601x1 Type=Int16, ColorInterp=Undefined NoData Value=-32768 Metadata: NETCDF_VARNAME=z
The grid is displayed correctly with GMT programs (grdimage, grdcontour), and is vertical-flipped with GDAL (GRASS, QGIS, gdal_translate).
I've attached a screenshot
Thanks José María
comment:3 by , 15 years ago
I've uploaded the screenshot twice. Please, someone, can delete the first copy, and this message? I can find the way to delete the image. Thanks, and sorry!
comment:4 by , 15 years ago
FWIW, the 'binary float' grid-registered version of the dataset loads into GRASS correctly using the r.in.bin module. see http://grass.osgeo.org/wiki/Global_datasets#ETOPO1
GDAL recognizes that binary float as "EHdr/ESRI .hdr Labelled" but incorrectly guesses Type=Byte, so r.in.gdal imports garbage.
Hamish
comment:5 by , 15 years ago
The EASE Grid projection (EPSG 3410) is also explained here: http://nsidc.org/data/atlas/epsg_3410.html (as PROJ4 parameters). Currently GDAL doesn't detect it in the netCDF files. Bruce pointed out that EASE Global undefined in the PROJ4 files:
# NSIDC EASE-Grid North <3408> +proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs <> # NSIDC EASE-Grid South <3409> +proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs <> # NSIDC EASE-Grid Global # Unable to translate coordinate system EPSG:3410 into PROJ.4 format.
The submitted file in this report should be theoretically detected as NSIDC EASE-Grid North.
comment:6 by , 15 years ago
I managed to extract metadata with Panoply netCDF viewer from http://www.giss.nasa.gov/tools/panoply/download_gen.html
java -jar Panoply.jar WARNING: Un-recognized grid_mapping : lambert_azimuthal_equal_area WARNING: Un-recognized grid_mapping : lambert_azimuthal_equal_area ... -------------- NetCDF file "NL.2004.nsidc0321v01.nc" netcdf file:/tmp/NL.2004.nsidc0321v01.nc { dimensions: time = UNLIMITED; // (46 currently) // (has coord.var) rows = 721; // (has coord.var) cols = 721; // (has coord.var) maxFilenameLength = 81; ncl4 = 10; variables: int date(time=46); :units = "yyyymmdd"; :long_name = "date"; char bpInfo(time=46, maxFilenameLength=81); :long_name = "SWE_file_used"; :_FillValue = ""; char scaInfo(time=46, maxFilenameLength=81); :long_name = "MODIS_percent_snow_covered_area_file_used"; :_FillValue = ""; short SWE(time=46, rows=721, cols=721); :standard_name = "lwe_thickness_of_surface_snow_amount"; :grid_mapping = "projection"; :ancillary_variables = "SCA bpInfo scaInfo"; :flag_meanings = " no_microwave_data grid_corners ocean permanent_ice_sheets_and_large_glaciers no_microwave_SWE_but_visible_SCA_gt_25_pcent"; :flag_values = -150s, -200s, -250s, -300s, -350s; // short :_FillValue = -200s; // short :information = " > 0 : SWE from deep microwave algorithm for 8-day period 0 : no snow -100 to -1 : -1 * SWE from shallow microwave algorithm -150 : missing microwave brightness temperatures -200 : fixed value for corners -250 : ocean -300 : permanent ice sheets and large glaciers -350 : no microwave SWE, but visible SCA > 25%"; :coordinates = "latitude longitude"; :units = "mm"; :long_name = "Snow Water Equivalent"; short SCA(time=46, rows=721, cols=721); :standard_name = "surface_snow_area_fraction"; :grid_mapping = "projection"; :ancillary_variables = "scaInfo"; :flag_meanings = " no_SCA_due_to_cloud_fill_night grid_corners ocean permanent_ice_sheets_and_large_glaciers"; :flag_values = -175s, -200s, -250s, -300s; // short :_FillValue = -200s; // short :information = " 0 : no snow 1 to 100 : percent MODIS snow-covered area -175 : no SCA, due to to cloud/fill/night -200 : fixed value for corners -250 : ocean -300 : permanent ice sheets and large glaciers"; :coordinates = "latitude longitude"; :units = "%"; :long_name = "MODIS snow-covered area (SCA)"; char projection(ncl4=10); :grid_id = "NL"; :reference = "EASE-Grid projection definition documention: http://nsidc.org/data/ease/ease_grid.html"; :NSIDC_mapx_grid_parameter_definition_file = "Nl.gpd"; :grid_mapping_name = "lambert_azimuthal_equal_area"; :longitude_of_projection_origin = 0.0f; // float :latitude_of_projection_origin = 90.0f; // float :false_easting = 0.0f; // float :false_northing = 0.0f; // float :projection_x_coordinate_origin = 360.0f; // float :projection_y_coordinate_origin = 360.0f; // float :scale_factor_at_projection_origin = 25.067526f; // float double time(time=46); :units = "days since -4713-01-01 12:00:00"; :long_name = "julian day"; :_CoordinateAxisType = "Time"; int rows(rows=721); :valid_range = 0, 720; // int :units = "count"; :standard_name = "projection_y_coordinate"; :long_name = "y-coordinate in EASE-Grid"; int cols(cols=721); :valid_range = 0, 720; // int :units = "count"; :standard_name = "projection_x_coordinate"; :long_name = "x-coordinate in EASE-Grid"; double latitude(rows=721, cols=721); :_FillValue = 1.0E20; // double :long_name = "latitude"; :units = "degrees_north"; :valid_range = -90.0f, 90.0f; // float :_CoordinateAxisType = "Lat"; double longitude(rows=721, cols=721); :_FillValue = 1.0E20; // double :long_name = "longitude"; :units = "degrees_east"; :valid_range = -180.0f, 180.0f; // float :_CoordinateAxisType = "Lon"; :creation_date = "Wed Feb 20 13:19:33 MST 2008"; :grid_id = "Nl"; :references = "Data set documentation: http://nsidc.org/data/nsidc-0321.html"; :comment = "Not set at this time"; :source = "Satellite observations from DMSP SSM/I and Terra MODIS"; :institution = "National Snow & Ice Data Center, Boulder, CO"; :title = "Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover"; :history = "Wed Feb 20 13:20:20 2008: ncrcat -O NL.2004001.nsidc0321v01.nc NL.2004009.nsidc0321v01.nc NL.2004017.nsidc0321v01.nc NL.2004025.nsidc0321v01.nc NL.2004033.nsidc0321v01.nc NL.2004041.nsidc0321v01.nc NL.2004049.nsidc0321v01.nc NL.2004057.nsidc0321v01.nc NL.2004065.nsidc0321v01.nc NL.2004073.nsidc0321v01.nc NL.2004081.nsidc0321v01.nc NL.2004089.nsidc0321v01.nc NL.2004097.nsidc0321v01.nc NL.2004105.nsidc0321v01.nc NL.2004113.nsidc0321v01.nc NL.2004121.nsidc0321v01.nc NL.2004129.nsidc0321v01.nc NL.2004137.nsidc0321v01.nc NL.2004145.nsidc0321v01.nc NL.2004153.nsidc0321v01.nc NL.2004161.nsidc0321v01.nc NL.2004169.nsidc0321v01.nc NL.2004177.nsidc0321v01.nc NL.2004185.nsidc0321v01.nc NL.2004193.nsidc0321v01.nc NL.2004201.nsidc0321v01.nc NL.2004209.nsidc0321v01.nc NL.2004217.nsidc0321v01.nc NL.2004225.nsidc0321v01.nc NL.2004233.nsidc0321v01.nc NL.2004241.nsidc0321v01.nc NL.2004249.nsidc0321v01.nc NL.2004257.nsidc0321v01.nc NL.2004265.nsidc0321v01.nc NL.2004273.nsidc0321v01.nc NL.2004281.nsidc0321v01.nc NL.2004289.nsidc0321v01.nc NL.2004297.nsidc0321v01.nc NL.2004305.nsidc0321v01.nc NL.2004313.nsidc0321v01.nc NL.2004321.nsidc0321v01.nc NL.2004329.nsidc0321v01.nc NL.2004337.nsidc0321v01.nc NL.2004345.nsidc0321v01.nc NL.2004353.nsidc0321v01.nc NL.2004361.nsidc0321v01.nc NL.2004.nsidc0321v01.nc"; :nco_openmp_thread_number = 1; // int }
Would be cool to read this through GDAL (hence, QGIS, GRASS etc...).
thanks, Markus
comment:7 by , 15 years ago
Frank,
I have created a new projection for LAEA and gdal seems to recognize it. Can someone test the new results?
Here is some CPLDebug output
OGRCT: Source: +proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs OGRCT: Target: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
Here is gdalinfo output:
gdalinfo NETCDF:"frmts/netcdf.old/NL.2004.nsidc0321v01.nc":SCA
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format Files: none associated Size is 721, 721 Coordinate System is: PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9108"]],
AXIS["Lat",NORTH], AXIS["Long",EAST], AUTHORITY["EPSG","4326"]],
PROJECTIONLambert_Azimuthal_Equal_Area, PARAMETER["latitude_of_center",90], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0]]
Origin = (-0.500000000000000,720.500000000000000) Pixel Size = (1.000000000000000,-1.000000000000000)
...
comment:8 by , 15 years ago
Denis,
it seems that the changes you've commited in r15689 and in r15690 have broken some autotests related to netcdf. See :
Running tests from gdrivers/netcdf.py TEST: netcdf_1 ... fail Checksum for band 1 in "NETCDF:"data/bug636.nc":tas" is 0, but expected 31621. TEST: netcdf_2 ... success TEST: netcdf_3 ... ERROR 1: netCDF scanline fetch failed: NetCDF: Numeric conversion not representable ERROR 1: IReadBlock failed at X offset 0, Y offset 0 ERROR 1: netCDF scanline fetch failed: NetCDF: Numeric conversion not representable ERROR 1: IReadBlock failed at X offset 0, Y offset 1 ERROR 1: netCDF scanline fetch failed: NetCDF: Numeric conversion not representable ... ERROR 1: netCDF scanline fetch failed: NetCDF: Numeric conversion not representable ERROR 1: IReadBlock failed at X offset 0, Y offset 100 fail Wrong min or max. TEST: netcdf_4 ... fail Checksum for band 3 in "NETCDF:data/foo_5dimensional.nc:temperature" is 613, but expected 1218.
follow-up: 10 comment:9 by , 15 years ago
Markus,
I have checked in an new version where I corrected the IReadBlock problem. I was able to display the data in OpenEV. I tried to do a gdalwarp on the file and it seems to work but the result did not satisfy me. I am not sure what to do at this point. Everything seem to be fine for your projection.
Give it a try.
wget ftp://sidads.colorado.edu/pub/DATASETS/snow/nsidc0321v01/north/NL.2004.nsidc0321v01.nc gdalinfo NETCDF:NL.2004.nsidc0321v01.nc:SCA Size is 721, 721 Coordinate System is: PROJCS["LAEA (WGS84) ", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9108"]], AXIS["Lat",NORTH], AXIS["Long",EAST], AUTHORITY["EPSG","4326"]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["latitude_of_center",90], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0]] ... OGRCT: Source: +proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs OGRCT: Target: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs Corner Coordinates: Upper Left ( -0.500, 720.500) (179d57'36.86"W, 89d59'36.78"N) Lower Left ( -0.5000000, -0.5000000) ( 45d 0'0.00"W, 89d59'59.98"N) Upper Right ( 720.500, 720.500) (135d 0'0.00"E, 89d59'27.16"N) Lower Right ( 720.500, -0.500) ( 89d57'36.86"E, 89d59'36.78"N) Center ( 360.000, 360.000) (135d 0'0.00"E, 89d59'43.59"N) ...
Denis
comment:10 by , 15 years ago
Denis,
I get the same result with gdalinfo which looks fine now. But:
Replying to dnadeau:
I tried to do a gdalwarp on the file and it seems to work but the result did not satisfy me.
what was exactly the command line here? I obtain only empty layers when using gdalwarp.
Markus
comment:11 by , 15 years ago
Denis,
the remaining problem might be that only row/col values are written instead of LAEA coodinates:
Corner Coordinates: Upper Left ( -0.500, 720.500) (179d57'36.86"W, 89d59'36.78"N) Lower Left ( -0.5000000, -0.5000000) ( 45d 0'0.00"W, 89d59'59.98"N) Upper Right ( 720.500, 720.500) (135d 0'0.00"E, 89d59'27.16"N) Lower Right ( 720.500, -0.500) ( 89d57'36.86"E, 89d59'36.78"N) Center ( 360.000, 360.000) (135d 0'0.00"E, 89d59'43.59"N)
If I take a country world map (FreeGIS world map) and reproject to EASE North, I obtain something different:
ogrinfo -so countries_simpl.shp countries_simpl | grep 'GEOGCS\|Extent' Extent: (-179.999900, -89.999900) - (179.999900, 83.609581) GEOGCS["wgs84", ogr2ogr -t_srs epsg:3408 world_ease_north.shp countries_simpl.shp ogrinfo -so world_ease_north.shp world_ease_north INFO: Open of `world_ease_north.shp' using driver `ESRI Shapefile' successful. Layer name: world_ease_north Geometry: Polygon Feature Count: 2461 Extent: (-12605642.317715, -12579528.665564) - (12482453.292396, 12742455.999976) Layer SRS WKT: PROJCS["NSIDC EASE-Grid North", GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere", DATUM["Mean_Sea_Level", SPHEROID["International_1924_Authalic_Sphere",6371228,0]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["latitude_of_center",90], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1]] cat: Real (11.0) vmap: String (80.0) fibs: String (80.0) iso3166: String (80.0) name: String (80.0)
The rendering in QGIS look ok. Could it be that the netCDF driver detects well the projection now but doesn't reproject the pixels yet?
Markus
comment:12 by , 12 years ago
The projection is in fact EPSG:3408
more info here
http://nsidc.org/data/ease/ease_grid.html#nh_map
comment:13 by , 9 years ago
By reading the comments the issue was due to difficulties in finding the projection and it was fixed some 6 years ago. I found that I was unable to confirm that. The data download links in this ticket are dead now and image formats used for delivering the snow and ice data are changed during these years. I found data from ftp://sidads.colorado.edu/pub/DATASETS/nsidc0046_weekly_snow_seaice/ but all the files are simple binary grids. Should they be wrapped somehow into NetCDF for making them readable for GDAL.
So unfortunately I can't say what is the status of this ticket. I hope that some of the cc:ed recipients can confirm that issue is fixed and this ticket can be closed.
comment:14 by , 6 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
Closing assuming everything is fixed. Reopen with more elements if not
EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover flipped east/west