gdal_translate returning PROJCS "unnamed" for lambert_cylindrical_equal_area

[brodzik]

I am using gdal_translate to convert netCDF data to geotiff for files in EASE-Grid 2.0 (equal area projections using WGS84 ellipsoid). I have successfully defined CF metadata for my Northern and Southern Azimuthal equal area files so that this command:

gdal_translate -of GTiff -b 1 NETCDF:"EASE2_N12.5km.F13_SSMI.1997061.37H.M.SIR.CSU.v0.1.nc":TB_num_samples N.tif

produces a geotiff that I can successfully read and reproject in ArcMap, for example. Also, gdalinfo returns

PROJCS["LAEA (WGS84) "

along with the other metadata I expect for GEOGCS and DATUM.

However, when I run the same gdal_translate command on my cylindrical equal-area file:

gdal_translate -of GTiff -b 1 NETCDF:"EASE2_T12.5km.F13_SSMI.1997061.85V.D.SIR.CSU.v0.1.nc":TB_num_samples T.tif

It returns PROJCS["unnamed"

but I would expect it to return "CEA (WGS84)" or something of the like.

I did notice an earlier ticket

​https://trac.osgeo.org/gdal/ticket/4068

that seemed partially related to this issue, for the original EASE-Grid CEA; it was also returning "unnamed" PROJCS, along with unnamed GEOGCS and DATUM. Theoretically, the EASE-Grid 2.0 definition I am using should have eliminated the projection spheroid issues, since I am using WGS84 for the projection.

I am using gdal 2.0.0 on OSX. I have put sample .nc files for North and cylindrical data, and the resulting gdal_translate *.tif files here:

​ftp://sidads.colorado.edu/pub/incoming/brodzik/cetb_samples/

Thank you for your attention.

[Even Rouault]

In 34881:

netCDF: add support for reading SRS from srid attribute when it exists and has content like urn:ogc:def:crs:EPSG::XXXX (fixes #6613)

[Even Rouault]

For EASE2_N12.5km.F13_SSMI.1997061.37H.M.SIR.CSU.v0.1.nc and LAEA, GDAL has a hard coded case to put "LAEA (WGS84)" for some reason, but generally there's no way to retrieve the projection citation from the CF attributes. I've added support for reading the SRS from the srid attribute that contains the EPSG code and enables GDAL to return a more expressive SRS.

[brodzik]

Evan,

I have been testing the gdal trunk for the changes you made, last time I pulled it was about a week ago. I am using ubuntu since it's easier to build from trunk this way. When I run gdal_translate, I see what looks like fine metadata, projection and geo CS look correct. However, ArcMap (10.4) is now throwing a nasty error when I try to read the geotiffs:

"The following data sources you added are missing spatial reference information. This data can be drawn in ArcMap but cannot be projected"

ArcMap does draw the data image, but sure enough the frame says Current Coordinate System is "No coordinate system".

I am putting examples of my .nc and the resulting geotiffs (from trunk gdal_translate) here:

​ftp://sidads.colorado.edu/pub/incoming/brodzik/testing_gdal/

and I am also trying to figure out how to submit an issue with ArcMap, but if you have any insight on what I might need to do to get ArcMap to understand the PROJCS information, I would be obliged.

Mary Jo

[Even Rouault]

Mari Jo,

my guess would be that ArcMap doesn't support the EPSG:6941 CRS encoded in the ProjectedCSTypeGeoKey geotiff key. Perhaps it was introduced in a too recent revision of the EPSG database not taken yet into account by ArcMap ?

What you could try is to encode the expanded definition of the CRS without the EPSG code. For example with:

gdal_translate NSIDC-0630-EASE2_N25km-F13_SSMI-1997061-19H-E-GRD-CSU-v1.0.nc.TB.tif out.tif -a_srs "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

[brodzik]

ESRI email regarding issue with spaces in PROJCS

[brodzik]

Frank,

Thanks for your message. I did the test you recommended using the proj4 string specifically, and it does produce a file that ArcMap opens and reprojects without error.

I have also been in touch with ESRI about the problem. They are suggesting that the problem is the spaces in the PROJCS name. They sent me the attached message which might make more sense to you than it does to me. They offered to speak by phone if that would be helpful and I agreed. We do not yet have a time arranged for the phone call.

When Evan addressed my last gdal issue by accessing the EPSG metadata directly I thought it was a great answer (since I went to all the trouble to get the EPSG definitions correct back in 2013). But now it sounds unfortunately like ESRI can't handle spaces in the ProjectedCRS name, is that how you read their message, also?

I will try to get a specific answer from ESRI when I speak to them. Would you consider taking part in that phone call?

Mary Jo

[Even Rouault]

Mary Jo, I think there's a confusion: the previous message was from me (Even), not Frank. I do not really believe in this theory of spaces (just my gut feeling. I can be wrong), as those texts are just informational and should normally not being used to get the projection. The point about .prj has little to do with GeoTIFF, especially since ESRI .prj does not include the EPSG code. Anyway, I've attached a crafted file (fake_6931_without_space.tif) where I manually replaced the spaces with underscore, so you can test it if it opens better in ArcMap (I think it will make no difference).

$ listgeo fake_6931_without_space.tif
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         440720            3751320           0                
      ModelPixelScaleTag (1,3):
         60                60                0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,35): "WGS_84___NSIDC_EASE_Grid_2_0_North"
      GeogCitationGeoKey (Ascii,7): "WGS_84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      ProjectedCSTypeGeoKey (Short,1): Unknown-6931
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

It should open centered at 173d17'21.40"E, 55d38'16.55"N. But my theory is that EPSG:6931 is too new for their database, rather than an issue with spaces. I don't have that much free time unfortunately to take part to a phone call.

[Even Rouault]

"gdal_translate byte.tif test.tif -a_srs EPSG:6931" and replacing spaces in geotiff keys by underscores in hexadecimal editor

[Even Rouault]

I dug into the GeoTIFF specification, and in ​http://web.archive.org/web/20160326103544/http://www.remotesensing.org/geotiff/spec/geotiff2.4.html#2.4, one can see the following example, which confirms that spaces are not forbidden, at least, by the spec :

  GeoAsciiParamsTag(34737)=("Custom File|My Geographic|")

[brodzik]

Even, Sorry about the confusion (and for misspelling your name)!

Thanks for the faked-out example. I am having trouble today connecting to the machine where we have ArcMap installed, I will test it Monday. In the meantime, I will pass it on to ESRI and proceed with the phone call to work with them. I really appreciate how responsive you are. I know my questions are complicated.

Mary Jo

[brodzik]

Even, I just got off the phone with ESRI tech support. We're trying to get to the bottom of this and she stepped me through some of the way that ArcMap is using the inputs. Turns out that they are using a very recent version of the EPSG data base and their defined projections do include the EPSG:6931-3 projections that I want. And we can see this in my version of ArcMap.

But, when ESRI imports the EPSG information, they have internal system restrictions on the projection name and they are doing some name-mangling to make their lookup work correctly. I want to help them run another test with the files that gdal is creating. Do you know of a way to change only the value of the string that's returned at the beginning of the PROJCS element in the gdalinfo call on the geotiffs? I'd like to produce another fake file that just has that string set to the value that ESRI mangled in their lookup, but with everything else you looked up with my EPSG srid. If it's something I can do with gdal tools and you can just point me to the right tool/options, that would be great.

If that's all they need and all the rest works as-is, I don't think I have any leverage to argue with ESRI that mangling's a bad idea. On a purely practical level, I may just have to fill out a change request with EPSG to change the ProjectedCRS name string to something they ESRI won't mangle. Bit of a moving target, but I'd like to give it a try.

Thanks, Mary Jo

[Even Rouault]

I fail to understand why ESRI doesn't just look at the content of the ProjectedCSTypeGeoKey key. This is where the EPSG code is. They shouldn't need anything else, and GDAL doesn't need anything else. How come do they expect interoperability with other software if they apply obscure rules ? Do they have an "official" document to explain their expectations (so to have possibly ways of exporting "ArcMap compatible" GeoTIFF) ?

As far as what you ask :

gdalsrsinfo -o wkt in.tif > wkt.txt
gdal_translate -a_srs wkt.txt in.tif out.tif

[brodzik]

Yes, I totally agree. Why make a "special" lookup table that no one else knows how to map to the authority? I will try to press this point with them. Yesterday I got the idea that the original tech support person is going to hand off my issue to someone closer to the technical side of things. Hopefully I can make the more general case with the new person, but it does feel awfully lonely trying to tell a big company like them how to do their job. Presumably the first step in influencing without authority is gaining some credibility with them. I don't really have anything to offer them in exchange (except for what I say to affect their reputation, which is pretty small coin).

Anyway, thanks for the gdalsrsinfo and gdal_translate recipe, that worked great. GDAL tools are wonderful.

I had to edit both the projection name from the original EPSG name and also remove the ",AUTHORITY["EPSG","6931"]" from the end of the WKT for it to work. Now when I run gdalinfo on the new tif file, the projection "looks" like their ESRI-mangled version.

And ArcMap opens and reprojects these ESRI-mangled .tif versions without throwing an error. I will take this back to them and try to argue against arbitrary, ESRI-defined lookup tables.

[Even Rouault]

Note that if you remove AUTHORITY["EPSG","6931"], then the geotiff keys are completely different than with it. They will be closer to what you get with -a_srs "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

You can see the effect with the listgeo utility

$ gdalsrsinfo -o wkt_simple EPSG:6931 > wkt_no_epsg.txt
$ gdal_translate in.tif out.tif -a_srs wkt_no_epsg.txt
$ listgeo out.tif
[...]
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,35): "WGS 84 / NSIDC EASE-Grid 2.0 North"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257223563    
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_LambertAzimEqualArea
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      ProjCenterLongGeoKey (Double,1): 0                
      ProjCenterLatGeoKey (Double,1): 90               
      End_Of_Keys.

whereas

$ gdal_translate in.tif out.tif -a_srs EPSG:6931
$ listgeo out.tif
[...]
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,35): "WGS 84 / NSIDC EASE-Grid 2.0 North"
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      ProjectedCSTypeGeoKey (Short,1): Unknown-6931
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.