Folks, 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.
Andreas, could you file a bug on this issue? Best regards, Frank On Wed, Apr 24, 2013 at 5:10 PM, David Shean <[email protected]> wrote: > 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 <[email protected]> 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 <[email protected]>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 >> [email protected] >> http://lists.osgeo.org/mailman/listinfo/gdal-dev >> > > > > -- > > ---------------------------------------+-------------------------------------- > I set the clouds in motion - turn up | Frank Warmerdam, > [email protected] > light and sound - activate the windows | http://pobox.com/~warmerdam > and watch the world go round - Rush | Geospatial Software Developer > _______________________________________________ > gdal-dev mailing list > [email protected] > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > > -- ---------------------------------------+-------------------------------------- I set the clouds in motion - turn up | Frank Warmerdam, [email protected] light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush | Geospatial Software Developer
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
