<div dir="ltr">Folks,<div><br></div><div style>It seems I neglected to review the output from the gdalinfo I posted before closely. I now notice it reports the latitude_of_origin as 0 isntead of -71 as it ought to for EPSG:3031. I believe this is a bug in the GeoTIFF code for reconsituting an SRS. On the other hand my "gdalsrsinfo EPSG:3031" does produce the latitude_of_origin as -71 as desired. </div>
<div style><br></div><div style>Andreas, could you file a bug on this issue? </div><div style><br></div><div style>Best regards,</div><div style>Frank</div><div style><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">
On Wed, Apr 24, 2013 at 5:10 PM, David Shean <span dir="ltr"><<a href="mailto:dshean@gmail.com" target="_blank">dshean@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word">Frank and Andreas,<br><div>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. <div>
<br></div><div>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.<div><br></div>
<div>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.</div>
<span class="HOEnZb"><font color="#888888"><div>-David</div></font></span><div><div class="h5"><div><br></div><div>On Apr 24, 2013, at 10:38 AM, Frank Warmerdam <<a href="mailto:warmerdam@pobox.com" target="_blank">warmerdam@pobox.com</a>> wrote:</div>
<div><div><br><blockquote type="cite"><div dir="ltr">Andreas,<div><br></div><div>I have downloaded test.3031.tif and run gdalinfo (trunk) against it and get the right results:</div><div><br></div><div><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 <a href="tel:%28-2054532.575" value="+12054532575" target="_blank">(-2054532.575</a>, 1166857.038) (119d35'38.75"W, 49d32' 9.40"N)</div>
<div>
Lower Left <a href="tel:%28-2054532.575" value="+12054532575" target="_blank">(-2054532.575</a>, 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 <a href="tel:%28-2029616.928" value="+12029616928" target="_blank">(-2029616.928</a>, 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>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><br></div><div>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><br></div><div>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><br></div><div><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" target="_blank">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>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><br></div><div><br></div><div>Best regards,<br></div><div>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" target="_blank">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>
_______________________________________________<br>gdal-dev mailing list<br><a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">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></blockquote>
</div><br></div></div></div></div></div></div></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>