[gdal-dev] gdal re-projeciton problem with Polar Stereographic projections
Frank Warmerdam
warmerdam at pobox.com
Wed Apr 24 17:41:20 PDT 2013
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 <dshean at gmail.com> 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 <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
>
>
>
--
---------------------------------------+--------------------------------------
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/143dac7c/attachment-0001.html>
More information about the gdal-dev
mailing list