Conversion to Proj for USGS

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Conversion to Proj for USGS

Ben Supnik
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.

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["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]]"
       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)

--
Scenery Home Page: http://scenery.x-plane.com/
Scenery blog: http://www.x-plane.com/blog/
Plugin SDK: http://www.xsquawkbox.net/xpsdk/
X-Plane Wiki: http://wiki.x-plane.com/
Scenery mailing list: [hidden email]
Developer mailing list: [hidden email]
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion to Proj for USGS

Even Rouault-2
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
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff