[Gdal-dev] eqc warping problem

[Frank Warmerdam]

On 10/2/05, Jianing Chen <jchenwiz at yahoo.com> wrote:
> Dear masters,
>
> I am running into a problem to warp images into eqc
> projection using lat_ts=37.0. I have attached the
> original with the warped. Please help out.
> Command line:
> gdalwarp.exe -t_srs "+proj=eqc +lat_ts=37.0 +lat_0=0.0
> +lon_0=0.0 +ellps=WGS84 +datum=WGS84 " -co
> "INTERLEAVE=PIXEL" -rb f:\testp.tif
> f:\output\testeqc.tif

JC,

To state the problem, it seems that the output file is created
with not quite accurate bounds.  It is essentially 20 meters south
of where it ought to be.   It helps to state why you think the result
is wrong!  Let me know if I guess wrong about the issue you are raising.

Looking at the data, you are reprojecting a point in Arizona State Plane
to Equidistant Conic.  Ostensibly you are transforming from NAD83 to
WGS84, but actually Equidistant Conic seems to be one of those
projections that only supports a spherical form.   If I try just going from
WGS84 to eqc coordinates with the PROJ.4 cs2cs command I see:

warmerda at gdal2200[92]% cs2cs +proj=latlong +datum=WGS84 +to +proj=eqc
+lat_ts=37.0 +lat_0=0.0 +lon_0=0.0 +ellps=WGS84 +datum=WGS84
-117 33
-10401732.70    3654001.32 -6305.98


Of particular note is the "Z" offset of over six kilometers. That is the
difference between the WGS84 ellipsoid and the sphere.

When gdalwarp works with coordinates it drops the Z.  This causes a
noticable difference between going "forward" from Arizona state plane
to eqc (at elevation 0) and going back from eqc to Arizona state plane
at elevation zero.

In short, reprojecting between a spherical projection like arizona
state plane and a spherical projection like eqc with gdalwarp since
gdalwarp ignores the substantial Z differences.

By the way, you can find out what projections PROJ.4 treats as
spherical only (vs. also supporting an ellipsoidal solution) using
the "proj -lP" command.

...
eqc : Equidistant Cylindrical (Plate Caree)
        Cyl, Sph
        lat_ts=
eqdc : Equidistant Conic
        Conic, Sph&Ell
        lat_1= lat_2=
...

Since eqc just has "Sph" and eqdc has "Sph&Ell" we see that
eqc is spherical only while eqdc also supports an ellipsoidal solution.

I am less sure what practical solution to offer.  I think what I may need to
do at some point is allow a mechanism such that PROJ.4's pj_transform()
method would not attempt to apply ellipsoid->sphere transformation via
geocentric coordinates.

Best regards,
--
---------------------------------------+--------------------------------------
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 Programmer for Rent