[gdal-dev] gdal re-projeciton problem with Polar Stereographic projections
David Shean
dshean at gmail.com
Wed Apr 24 17:10:36 PDT 2013
Frank and Andreas,
I've experienced similar issues for years when using EPSG codes for polar stereographic projections (e.g., EPSG:3031, EPSG:3413). Thanks for reporting Andreas.
My workaround for this issue is to provide the full proj4 string to gdalwarp instead of the EPSG code. This string can be extracted with gdalsrsinfo utility or the osr ImportFromEPSG method.
I just performed a quick test with trunk gdalwarp (GDAL 1.10.0, released 2013/04/13) with -t_srs EPSG:3413, and experienced the same issues described by Andreas. Using a UTM code for -t_srs EPSG:32622 works fine, with appropriate PROJCS tags in the output gtiff.
-David
On Apr 24, 2013, at 10:38 AM, Frank Warmerdam <warmerdam at pobox.com> wrote:
> 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
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130424/61133c8f/attachment.html>
More information about the gdal-dev
mailing list