6

I've been trying to use GDAL to convert from an arbitrary projection/datum into Equirectangular (aka Plate Caree, aka Equidistant Cylindrical, aka Geographic, aka Lat/Long) with WGS84 datum. All of my input files have an associated PRJ file that I've been using with GDAL with this statement:

gdalwarp -s_srs Input.prj -t_srs Output.prj Input.xyz Output.tif

So far, I've had mixed results. I've been able to successfully transform a Lambert Equal Area projection in WGS84 to Geographic/WGS84, but performing a transform on LEA in RT90 projection causes an error of ~1km.

Stranger still, I tried to transform two Sinusoidal projections in RT90 and in WGS84, both of which resulted in an error of 16km! I suspect that these errors are largely due to some sort of ellipsoid/sphere discrepancy.

I'm very confused as to why these errors are so random, but more concerned that they are happening at all. Shouldn't gdalwarp take any ellipsoid information about the source and target from the PRJ files and transform between them?

What can I do to correct my issue?

|improve this question
  • please provide the contents of the Input.prj and the Output.prj for the various cases of success and failure that you report – mdsumner Dec 1 '10 at 20:12
  • 1
    This is a cross posting: lists.osgeo.org/pipermail/gdal-dev/2010-November/026929.html – mdsumner Dec 1 '10 at 20:25
  • Its important to remember that resampling and reprojection of rasters is a lossy process, even after you sort out your projection issues. – scw Dec 1 '10 at 21:33
4

I've solved my problem.

RT90 is not a defined datum in Proj4. As such, GDAL just up and ignores the datum entry and applies the ellipsoid and projection correction.

It'd be nice if a warning was thrown at some point ("WARNING: datum 'RT90' not recognized" would be enough, but I digress).

The sinusoidal projection in WGS84 datum was transforming incorrectly because the Sinusoidal projection wasn't being set properly when export from Global Mapper. When I hand-rolled my own SRS to be "+proj=sinu +datum=WGS84" everything worked out fine.

To get proj4 to list various types of datums/projections/ellipsoids that it supports:

proj -ld    //lists known datums
proj -le    //lists known ellipsoids
proj -l     //lists known projections (use -lP for an expanded view)

If you want GDAL to work with unknown datum types, you must specify "+towgs84=" in your source SRS, and you must also specify "+towgs84=" in your target SRS if your target datum is not WGS84. These tags must be followed by 3 or 7 parameter transformation values, which can be found with some google-fu. More info on this tag can be found here:

http://trac.osgeo.org/proj/wiki/GenParms#towgs84-DatumtransformationtoWGS84

Hopefully anyone that runs into a similar problem will now be able to do so.

|improve this answer
2

There are many cases where transformation between datums requires more complete information like the parameters for towgs84 - both for input and output - presumably you have left this out for the transformations that are failing, or you have somehow specified them incorrectly.

|improve this answer
  • I guess I'm just confused as to why GDAL doesn't have these values baked in. As per your earlier request, here are the Input and Output files: – Scott S Dec 2 '10 at 15:17
  • Input.prj (Sinusoidal/RT90): PROJCS["Sinusoidal",GEOGCS["Geographic Coordinate System",DATUM["RT90",SPHEROID["Bessel 1841",6377397.155,299.1528131060786]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]] – Scott S Dec 2 '10 at 15:18
  • Output.prj (Lat/Long, WGS84): GEOGCS["Geographic Coordinate System",DATUM["WGS84",SPHEROID["WGS84",6378137,298.257223560493]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] – Scott S Dec 2 '10 at 15:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy