<div dir="ltr">Andreas,<div><br></div><div style>I have downloaded test.3031.tif and run gdalinfo (trunk) against it and get the right results:</div><div style><br></div><div style><div>Driver: GTiff/GeoTIFF</div><div>Files: test.3031.tif</div>
<div>Size is 398, 398</div><div>Coordinate System is:</div><div>PROJCS["WGS 84 / Antarctic Polar Stereographic",</div><div> GEOGCS["WGS 84",</div><div> DATUM["WGS_1984",</div><div> SPHEROID["WGS 84",6378137,298.257223563,</div>
<div> AUTHORITY["EPSG","7030"]],</div><div> AUTHORITY["EPSG","6326"]],</div><div> PRIMEM["Greenwich",0],</div><div> UNIT["degree",0.0174532925199433],</div>
<div> AUTHORITY["EPSG","4326"]],</div><div> PROJECTION["Polar_Stereographic"],</div><div> PARAMETER["latitude_of_origin",0],</div><div> PARAMETER["central_meridian",0],</div>
<div> PARAMETER["scale_factor",1],</div><div> PARAMETER["false_easting",0],</div><div> PARAMETER["false_northing",0],</div><div> UNIT["metre",1,</div><div> AUTHORITY["EPSG","9001"]],</div>
<div> AUTHORITY["EPSG","3031"]]</div><div>Origin = (-2054532.574583648238331,1166857.038436320610344)</div><div>Pixel Size = (125.204255282867464,-125.204255282867464)</div><div>Metadata:</div><div>
AREA_OR_POINT=Area</div><div>Image Structure Metadata:</div><div> INTERLEAVE=BAND</div><div>Corner Coordinates:</div><div>Upper Left (-2054532.575, 1166857.038) (119d35'38.75"W, 49d32' 9.40"N)</div><div>
Lower Left (-2054532.575, 1117025.745) (118d31'56.70"W, 49d55' 6.39"N)</div><div>Upper Right (-2004701.281, 1166857.038) (120d12' 7.15"W, 50d13' 9.00"N)</div><div>Lower Right (-2004701.281, 1117025.745) (119d 7'36.07"W, 50d36'37.86"N)</div>
<div>Center (-2029616.928, 1141941.392) (119d21'49.67"W, 50d 4'21.54"N)</div><div>Band 1 Block=398x10 Type=UInt16, ColorInterp=Gray</div><div><br></div><div><br></div><div style>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). </div>
<div style><br></div><div style>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. </div>
<div style><br></div><div style>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. </div>
<div style><br></div><div style><br></div></div><div class="gmail_extra"><div class="gmail_extra">==== //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 ====</div>
<div class="gmail_extra">137,138c137,146</div><div class="gmail_extra">< // Polar A case</div><div class="gmail_extra">< latitudeTrueScale = (latitudeOfOrigin < 0) ? -PI/2 : PI/2;</div><div class="gmail_extra">
---</div><div class="gmail_extra">> // If we do not have standard_parallel, but we do have a reasonable</div><div class="gmail_extra">> // value for latitude of origin, then use that as the latitude of</div>
<div class="gmail_extra">> // true scale as documented at:</div><div class="gmail_extra">> // <a href="http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html">http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html</a></div>
<div class="gmail_extra">> if (abs(latitudeOfOrigin) > PI / 4) {</div><div class="gmail_extra">> latitudeTrueScale = latitudeOfOrigin;</div><div class="gmail_extra">> } else {</div>
<div class="gmail_extra">> // Polar A case</div><div class="gmail_extra">> latitudeTrueScale = (latitudeOfOrigin < 0) ? -PI / 2 : PI / 2;</div><div class="gmail_extra">> }</div>
<div class="gmail_extra">155d162</div><div class="gmail_extra">< this.latitudeOfOrigin = (southPole) ? -(PI/2) : +(PI/2);</div><div class="gmail_extra">161c168</div><div class="gmail_extra">< if (abs(latitudeTrueScale - PI/2) >= EPSILON) {</div>
<div class="gmail_extra">---</div><div class="gmail_extra">> if (abs(latitudeTrueScale - PI / 2) >= EPSILON) {</div><div class="gmail_extra">386c393,397</div><div class="gmail_extra">< y = latitudeOfOrigin;</div>
<div class="gmail_extra">---</div><div class="gmail_extra">> if (southPole) {</div><div class="gmail_extra">> y = -(PI / 2);</div><div class="gmail_extra">> } else {</div>
<div class="gmail_extra">> y = (PI / 2);</div><div class="gmail_extra">> }</div><div><br></div>You can also manually manipulate the WKT coordinate system definition to get it to work with GeoTools. As a hackaround I used:</div>
<div class="gmail_extra"><div class="gmail_extra"><br></div><div class="gmail_extra">> if (wkt.contains("Polar_Stereographic")) {</div><div class="gmail_extra">> wkt = wkt</div><div class="gmail_extra">
> .replace("latitude_of_origin", "Standard_Parallel_1")</div><div class="gmail_extra">> .replace("Polar_Stereographic", "Polar Stereographic (variant B)")</div>
<div class="gmail_extra">> .replace(",PARAMETER[\"scale_factor\",1]", "");</div><div class="gmail_extra">> }</div><div><br></div><div style>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. </div>
<div style><br></div><div style><br></div><div style>Best regards,<br></div><div style>Frank</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Apr 24, 2013 at 9:54 AM, Cziferszky, Andreas <span dir="ltr"><<a href="mailto:ancz@bas.ac.uk" target="_blank">ancz@bas.ac.uk</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hi list,<br>
<br>
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<br>
<br>
gdalwarp -s_srs EPSG:3395 -t_srs EPSG:3031 test.3395.tif test.3031.tif<br>
<br>
generates a GeoTIFF whose map projection info cannot be read using ArcMap (v10.1 in our case). For the out image gdalinfo is reporting<br>
<br>
>gdalinfo test.3031.tif<br>
Driver: GTiff/GeoTIFF<br>
Files: test.3031.tif<br>
Size is 398, 398<br>
Coordinate System is:<br>
LOCAL_CS["WGS 84 / Antarctic Polar Stereographic",<br>
GEOGCS["WGS 84",<br>
DATUM["WGS_1984",<br>
SPHEROID["WGS 84",6378137,298.257223563,<br>
AUTHORITY["EPSG","7030"]],<br>
AUTHORITY["EPSG","6326"]],<br>
PRIMEM["Greenwich",0],<br>
UNIT["degree",0.0174532925199433],<br>
AUTHORITY["EPSG","4326"]],<br>
AUTHORITY["EPSG","3031"],<br>
UNIT["metre",1]]<br>
....<br>
<br>
Two items seem to be wrong/missing:<br>
1) The LOCAL_CS definition instead of an expected PROJCS definition<br>
2) No projection parameters like PROJECTION["Polar Stereographic..."], PARAMETER["central_meridian",0], UNIT["metre",1,..].. exist, only the AUTHORITY including EPSG code 3031.<br>
<br>
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.<br>
<br>
Files are available at <a href="ftp://ftp.bas.ac.uk/ancz/test.3395.tif" target="_blank">ftp://ftp.bas.ac.uk/ancz/test.3395.tif</a> and <a href="ftp://ftp.bas.ac.uk/ancz/test.3031.tif" target="_blank">ftp://ftp.bas.ac.uk/ancz/test.3031.tif</a><br>
<br>
Any help appreciated.<br>
Andreas<br>
<br>
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.<br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>---------------------------------------+--------------------------------------<br>I set the clouds in motion - turn up | Frank Warmerdam, <a href="mailto:warmerdam@pobox.com" target="_blank">warmerdam@pobox.com</a><br>
light and sound - activate the windows | <a href="http://pobox.com/~warmerdam" target="_blank">http://pobox.com/~warmerdam</a><br>and watch the world go round - Rush | Geospatial Software Developer<br>
</div></div>