[gdal-dev] gdal re-projeciton problem with Polar Stereographic projections
Frank Warmerdam
warmerdam at pobox.com
Wed Apr 24 10:38:57 PDT 2013
Andreas,
I have downloaded test.3031.tif and run gdalinfo (trunk) against it and get
the right results:
Driver: GTiff/GeoTIFF
Files: test.3031.tif
Size is 398, 398
Coordinate System is:
PROJCS["WGS 84 / Antarctic Polar Stereographic",
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"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3031"]]
Origin = (-2054532.574583648238331,1166857.038436320610344)
Pixel Size = (125.204255282867464,-125.204255282867464)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-2054532.575, 1166857.038) (119d35'38.75"W, 49d32' 9.40"N)
Lower Left (-2054532.575, 1117025.745) (118d31'56.70"W, 49d55' 6.39"N)
Upper Right (-2004701.281, 1166857.038) (120d12' 7.15"W, 50d13' 9.00"N)
Lower Right (-2004701.281, 1117025.745) (119d 7'36.07"W, 50d36'37.86"N)
Center (-2029616.928, 1141941.392) (119d21'49.67"W, 50d 4'21.54"N)
Band 1 Block=398x10 Type=UInt16, ColorInterp=Gray
I don't know if this was something fixed after 1.9.2 or if the problem is
the configuration of your coordinate system support files (ie pcs.csv).
I'd encourage you to upgrade. to 1.10.0 (the RC3 is available).
With regard to using this in GeoServer, somewhat strangely I am in the
midst of preparing a patch for GeoTools so that the WKT representation of
this produced by GDAL will work with GeoTools (and so GeoServer). I
haven't yet filed the bug against GeoTools for this since I am still
working to get this through review internally here.
The gist of the change is that GeoTools does not recognise the
latitude_of_origin as also being the latitude of true scale - it uses the
pole for that. Here is a patch against (I think) GeoTools 8.1.
====
//depot/google3/third_party/java_src/geotools/gt_referencing/java/org/geotools/referencing/operation/projection/PolarStereographic.java#1
-
/google/src/cloud/warmerdam/ee2/google3/third_party/java_src/geotools/gt_referencing/java/org/geotools/referencing/operation/projection/PolarStereographic.java
====
137,138c137,146
< // Polar A case
< latitudeTrueScale = (latitudeOfOrigin < 0) ? -PI/2 : PI/2;
---
> // If we do not have standard_parallel, but we do have a
reasonable
> // value for latitude of origin, then use that as the
latitude of
> // true scale as documented at:
> //
http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html
> if (abs(latitudeOfOrigin) > PI / 4) {
> latitudeTrueScale = latitudeOfOrigin;
> } else {
> // Polar A case
> latitudeTrueScale = (latitudeOfOrigin < 0) ? -PI / 2 : PI /
2;
> }
155d162
< this.latitudeOfOrigin = (southPole) ? -(PI/2) : +(PI/2);
161c168
< if (abs(latitudeTrueScale - PI/2) >= EPSILON) {
---
> if (abs(latitudeTrueScale - PI / 2) >= EPSILON) {
386c393,397
< y = latitudeOfOrigin;
---
> if (southPole) {
> y = -(PI / 2);
> } else {
> y = (PI / 2);
> }
You can also manually manipulate the WKT coordinate system definition to
get it to work with GeoTools. As a hackaround I used:
> if (wkt.contains("Polar_Stereographic")) {
> wkt = wkt
> .replace("latitude_of_origin", "Standard_Parallel_1")
> .replace("Polar_Stereographic", "Polar Stereographic (variant
B)")
> .replace(",PARAMETER[\"scale_factor\",1]", "");
> }
Lastly I will note that ESRI represents Polar Stereographic slightly
differently than GDAL. It uses standard_parallel_1 instead of
latitude_of_origin for the latitude of true scale. GDAL normally knows how
to make this adjustment when converting.
Best regards,
Frank
On Wed, Apr 24, 2013 at 9:54 AM, Cziferszky, Andreas <ancz at bas.ac.uk> wrote:
> Hi list,
>
> We are experiencing problems with GeoTIFF images generated by gdalwarp
> (v1.9.2) while map projecting towards a Polar Stereographic projection,
> particularly EPSG:3031 (although other PS map projections showed the same
> effect). The source image is a World Mercator (EPSG:3395) projected GeoTIFF
> which displays and reports map proejction info correctly in QGIS, ArcMap,
> and others. Running
>
> gdalwarp -s_srs EPSG:3395 -t_srs EPSG:3031 test.3395.tif test.3031.tif
>
> generates a GeoTIFF whose map projection info cannot be read using ArcMap
> (v10.1 in our case). For the out image gdalinfo is reporting
>
> >gdalinfo test.3031.tif
> Driver: GTiff/GeoTIFF
> Files: test.3031.tif
> Size is 398, 398
> Coordinate System is:
> LOCAL_CS["WGS 84 / Antarctic Polar Stereographic",
> 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"]],
> AUTHORITY["EPSG","3031"],
> UNIT["metre",1]]
> ....
>
> Two items seem to be wrong/missing:
> 1) The LOCAL_CS definition instead of an expected PROJCS definition
> 2) No projection parameters like PROJECTION["Polar Stereographic..."],
> PARAMETER["central_meridian",0], UNIT["metre",1,..].. exist, only the
> AUTHORITY including EPSG code 3031.
>
> In addition, it appears that some (not all) GeoTIFF files projected using
> gdalwarp (1.9.1) into Polar Stereographic (EPSG:3031) as above cause
> problems when publishing via the Geoserver admin GUI (v2.2.x). The problem
> appears when trying to save the raster layer's properties - the native SRS
> appears undefined, presumably because Geotools has had issues with the
> LOCAL_CS header, so that only "reproject native to declared" actually
> works. The other two options ("keep native" and "force declared") cause a
> stack trace, and failure to save the layer properties.
>
> Files are available at ftp://ftp.bas.ac.uk/ancz/test.3395.tif and
> ftp://ftp.bas.ac.uk/ancz/test.3031.tif
>
> Any help appreciated.
> Andreas
>
> This message (and any attachments) is for the recipient only. NERC is
> subject to the Freedom of Information Act 2000 and the contents of this
> email and any reply you make may be disclosed by NERC unless it is exempt
> from release under the Act. Any material supplied to NERC may be stored in
> an electronic records management system.
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam,
warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Software Developer
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130424/2a00beef/attachment-0001.html>
More information about the gdal-dev
mailing list