[Geotiff] Conversion to Proj for USGS

Sun May 8 05:50:35 EST 2016

Le mardi 10 mars 2015 00:51:38, Ben Supnik a écrit :
> Hi Y'all,
> 
> In testing the latest version of our app, one of our testers grabbed a
> NAIP GeoJP2K from the USGS, and found that it won't unproject.  (This
> isn't off topic, really! :-)
> 
> The spacial data for the GeoJP2K comes from a stub GeoTIFF file embedded
> inside, with an unused 1x1 pixel and GeoTIFF tags that set up the
> coordinate systems, etc.
> 
> Apparently the USGS is using "WGS_1984_Web_Mercator_Auxiliary_Sphere",
> which is apparently a cheap spherical mercator projector designed to
> match Google/Bing and other wide-scale commercial mapping services.
> 
> I looked at the code, and the problem appears to be that while there is
> a proj string that will unproject this, libgeotiff doesn't know how to
> extract it from the given geo tags.  I'm not sure I blame it - I see
> text-based identification of the projection and an ESRI PE string, but
> no EPSG tag or enumerated description of the projection.
> 
> My stupid questions:
> 
> - Is there a way to mechanically build the projection string from
> anything in this set of GeoTIFF tags?
> - If not, does anyone have any suggestions?  Our user base are not core
> GIS users, so having them hand-specify spatial parameters or edit the
> spatial information isn't really on the table.

Ben,

Yes, support for EPSG:3857 / Google WebMercator is a mess for 
interoperability, especially in GeoTIFF since it encoding it is not 
standardized, so different vendors do different things. Some will use a 
ProjectedCSTypeGeoKey (Short,1): Unknown-3857 (GDAL), some will use what you 
got (ESRI). 

In recent version of GDAL, I've added logic to understand some of the 
formulations, and in particular the one you mention, so the correct proj.4 
string "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 
+y_0=0 +k=1.0 +units=m +nadgrids=@null  +no_defs" will now be returned.  
geotiff_proj4.c from libgeotiff could potentially be tweaked also to do similar 
things (GDAL doesn't use it, but directly deal with the geo_normalize.h 
interface)

$ gdalinfo out.tif
Driver: GTiff/GeoTIFF
Files: out.tif
Size is 20, 20
Coordinate System is:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mercator_Auxiliary_Sphere"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    PARAMETER["Standard_Parallel_1",0.0],
    PARAMETER["Auxiliary_Sphere_Type",0.0],
    UNIT["Meter",1.0],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 
+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"]]
Origin = (-12481700.191124100238085,4300624.213151689618826)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-12481700.191, 4300624.213) (112d 7'30.07"W, 36d 0' 0.07"N)
Lower Left  (-12481700.191, 4300604.213) (112d 7'30.07"W, 35d59'59.55"N)
Upper Right (-12481680.191, 4300624.213) (112d 7'29.43"W, 36d 0' 0.07"N)
Lower Right (-12481680.191, 4300604.213) (112d 7'29.43"W, 35d59'59.55"N)
Center      (-12481690.191, 4300614.213) (112d 7'29.75"W, 35d59'59.81"N)
Band 1 Block=20x20 Type=Byte, ColorInterp=Gray

Best regards,

Even

> 
> The listgeo dump is listed below.
> 
> Thanks!
> Ben
> 
> Geotiff_Information:
>     Version: 1
>     Key_Revision: 1.0
>     Tagged_Information:
>        ModelTiepointTag (2,3):
>           0                 0                 0
>           -12481700.1911241 4300624.21315169  0
>        ModelPixelScaleTag (1,3):
>           1                 1                 0
>        End_Of_Tags.
>     Keyed_Information:
>        GTModelTypeGeoKey (Short,1): User-Defined
>        GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
>        GTCitationGeoKey (Ascii,50): "PCS Name =
> WGS_1984_Web_Mercator_Auxiliary_Sphere"
>        GeographicTypeGeoKey (Short,1): GCS_WGS_84
>        GeogCitationGeoKey (Ascii,13): "GCS_WGS_1984"
>        GeogAngularUnitsGeoKey (Short,1): Angular_Degree
>        GeogSemiMajorAxisGeoKey (Double,1): 6378137
>        GeogInvFlatteningGeoKey (Double,1): 298.257223563
>        PCSCitationGeoKey (Ascii,445): "ESRI PE String =
> PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM
> ["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwi
> ch",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliar
> y_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],P
> ARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAM
> ETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]" ProjLinearUnitsGeoKey
> (Short,1): Linear_Meter
>        End_Of_Keys.
>     End_Of_Geotiff.
> 
> GCS: 4326/WGS 84
> Datum: 6326/World Geodetic System 1984
> Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
> Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
> Projection Linear Units: 9001/metre (1.000000m)
> 
> Corner Coordinates:
> Upper Left    (-12481700.191, 4300624.213)
> Lower Left    (-12481700.191, 4300623.213)
> Upper Right   (-12481699.191, 4300624.213)
> Lower Right   (-12481699.191, 4300623.213)
> Center        (-12481699.691, 4300623.713)

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com