gdalwarp EPSG:32662 problem

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

gdalwarp EPSG:32662 problem

Stephen Woodbridge
Hi,

Sorry for cross posting. This seems like a proj4 issue but might be gdal
related.

I've run into a problem that I think is related to a change in the EPSG
definition for EPSG:32662

The symptom is that the image after gdalwarp reprojection is about 20km
south of where it should be.

old system is correctly aligned:

GDAL 1.4.2.0, released 2007/06/27
proj Rel. 4.5.0, 22 Oct 2006

# WGS 84 / Plate Carree
<32662> +proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +units=m +no_defs  <>

new system is shifted about 20km south:

GDAL 1.9.2, released 2012/10/08
proj Rel. 4.8.0, 6 March 2012

# WGS 84 / Plate Carree (deprecated)
<32662> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs  <>


http://spatialreference.org/ref/epsg/32662/
This indicates that the definition in proj 4.8.0 is wrong. Would the
+ellps=WGS84 account for this shift?

Any thoughts on this?

Thanks,
   -Steve
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Fwd: [gdal-dev] gdalwarp EPSG:32662 problem

Antti Castrén-2
Hello Steve,

You asked for any thoughts, so here are some.

The 20km shift could easily be caused by using spherical version in
one implementation of equirectangular projection and ellipsoidal
version in the other. Amount of shift depends on the latitude.

You can find some information on these pages:

http://www.remotesensing.org/geotiff/proj_list/equirectangular.html

http://comments.gmane.org/gmane.comp.gis.proj-4.devel/981

Read especially the comments of the latter.

Hope these help,

  Antti

2014-02-08 8:11 GMT+02:00 Stephen Woodbridge <[hidden email]>:

> Hi,
>
> Sorry for cross posting. This seems like a proj4 issue but might be gdal
> related.
>
> I've run into a problem that I think is related to a change in the EPSG
> definition for EPSG:32662
>
> The symptom is that the image after gdalwarp reprojection is about 20km
> south of where it should be.
>
> old system is correctly aligned:
>
> GDAL 1.4.2.0, released 2007/06/27
> proj Rel. 4.5.0, 22 Oct 2006
>
> # WGS 84 / Plate Carree
> <32662> +proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs  <>
>
> new system is shifted about 20km south:
>
> GDAL 1.9.2, released 2012/10/08
> proj Rel. 4.8.0, 6 March 2012
>
> # WGS 84 / Plate Carree (deprecated)
> <32662> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs  <>
>
>
> http://spatialreference.org/ref/epsg/32662/
> This indicates that the definition in proj 4.8.0 is wrong. Would the
> +ellps=WGS84 account for this shift?
>
> Any thoughts on this?
>
> Thanks,
>   -Steve
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: gdalwarp EPSG:32662 problem

support.mn
In reply to this post by Stephen Woodbridge
Hello,

the new version of proj.4 does not make any datum or ellipsoidal
shift if not both source and destination datums are defined.. the older
version did some other assumptions.. this is the most often reason
for some shifts between versions.

see also:
http://comments.gmane.org/gmane.comp.gis.proj-4.devel/6045

Regards: Janne.

---------------------------------------------------------------------------------------


Stephen Woodbridge [[hidden email]] kirjoitti:

> Hi,
>
> Sorry for cross posting. This seems like a proj4 issue but might be gdal
> related.
>
> I've run into a problem that I think is related to a change in the EPSG
> definition for EPSG:32662
>
> The symptom is that the image after gdalwarp reprojection is about 20km
> south of where it should be.
>
> old system is correctly aligned:
>
> GDAL 1.4.2.0, released 2007/06/27
> proj Rel. 4.5.0, 22 Oct 2006
>
> # WGS 84 / Plate Carree
> <32662> +proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
> +datum=WGS84 +units=m +no_defs  <>
>
> new system is shifted about 20km south:
>
> GDAL 1.9.2, released 2012/10/08
> proj Rel. 4.8.0, 6 March 2012
>
> # WGS 84 / Plate Carree (deprecated)
> <32662> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs  <>
>
>
> http://spatialreference.org/ref/epsg/32662/
> This indicates that the definition in proj 4.8.0 is wrong. Would the
> +ellps=WGS84 account for this shift?
>
> Any thoughts on this?
>
> Thanks,
>    -Steve
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: gdalwarp EPSG:32662 problem

Stephen Woodbridge
Hi Janne, Antti,

I gues a better question is using gdalwarp, what is the correct way to
project this image so it aligns. given the behavior has changed between
the two versions of proj4 that I listed below.

Reading the links people have posted explains the difference, but it is
not clear to me if changing the epsg file or the command line and how is
the right thing to do here.

Thanks,
   -Steve

My gdalwarp command is:

gdalwarp -rcs -ts 8800 6600 -s_srs EPSG:32662 -t_srs EPSG:4326
$temptiff1 $target

The gdalinfo on the image is:

root@u17269306:/maps/images.tmp# gdalinfo
MODIS.2014038.193758.fullpass.seadas_sst.hdf
Driver: HDF4/Hierarchical Data Format Release 4
Files: MODIS.2014038.193758.fullpass.seadas_sst.hdf
Size is 512, 512
Coordinate System is `'
Metadata:
   Collection Location=USF Institute for Marine Remote Sensing
   Data Source Type=Direct Broadcast
   Input File(s)=/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038193749.GEO,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038193749.L2_LAC_SST,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194019.GEO,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194019.L2_LAC_SST,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194249.GEO,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194249.L2_LAC_SST,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194519.GEO,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194519.L2_LAC_SST,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194749.GEO,
/tmp/s4p_run_map_modis_2.3_WSthzYwz/A2014038194749.L2_LAC_SST
   L3 Software Name=/usr/bin/map_modis.pl
   L3 Software Version=4.2
   Map Limit=-5, -125, 55, -45
   Map Projection Category=Proj 4.4.7
   Map Projection Name=eqc
   Processing Status=intermediate
   Proj Projection Parameters=#Equidistant Cylindrical (Plate Caree)
#       Cyl, Sph
#       lat_ts=
# +proj=eqc +ellps=WGS84
   Satellite=aqua
   Spatial Resolution=1km
   Temporal Resolution=pass
Subdatasets:
 
SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"MODIS.2014038.193758.fullpass.seadas_sst.hdf":0
   SUBDATASET_1_DESC=[6600x8800] seadas_sst (16-bit integer)
 
SUBDATASET_2_NAME=HDF4_SDS:UNKNOWN:"MODIS.2014038.193758.fullpass.seadas_sst.hdf":1
   SUBDATASET_2_DESC=[6600x8800] l2_flags (32-bit integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

On 2/8/2014 8:23 AM, [hidden email] wrote:

> Hello,
>
> the new version of proj.4 does not make any datum or ellipsoidal
> shift if not both source and destination datums are defined.. the older
> version did some other assumptions.. this is the most often reason
> for some shifts between versions.
>
> see also:
> http://comments.gmane.org/gmane.comp.gis.proj-4.devel/6045
>
> Regards: Janne.
>
> ---------------------------------------------------------------------------------------
>
>
> Stephen Woodbridge [[hidden email]] kirjoitti:
>> Hi,
>>
>> Sorry for cross posting. This seems like a proj4 issue but might be gdal
>> related.
>>
>> I've run into a problem that I think is related to a change in the EPSG
>> definition for EPSG:32662
>>
>> The symptom is that the image after gdalwarp reprojection is about 20km
>> south of where it should be.
>>
>> old system is correctly aligned:
>>
>> GDAL 1.4.2.0, released 2007/06/27
>> proj Rel. 4.5.0, 22 Oct 2006
>>
>> # WGS 84 / Plate Carree
>> <32662> +proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
>> +datum=WGS84 +units=m +no_defs  <>
>>
>> new system is shifted about 20km south:
>>
>> GDAL 1.9.2, released 2012/10/08
>> proj Rel. 4.8.0, 6 March 2012
>>
>> # WGS 84 / Plate Carree (deprecated)
>> <32662> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
>> +units=m +no_defs  <>
>>
>>
>> http://spatialreference.org/ref/epsg/32662/
>> This indicates that the definition in proj 4.8.0 is wrong. Would the
>> +ellps=WGS84 account for this shift?
>>
>> Any thoughts on this?
>>
>> Thanks,
>>     -Steve
>> _______________________________________________
>> Proj mailing list
>> [hidden email]
>> http://lists.maptools.org/mailman/listinfo/proj
>>
>
>
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: [gdal-dev] gdalwarp EPSG:32662 problem

Antti Castrén-2
In reply to this post by Stephen Woodbridge
2014-02-09 8:40 GMT+02:00 Andre Joost <[hidden email]>:

> Am 09.02.2014 00:42, schrieb Even Rouault:
>>
>> But Antti guess seems right. Instead of +ellps=WGS84 (or +datum=WGS84), if
>> you
>> play with the a (semi-major axis) and b (semi-minor axis) parameters, you
>> can
>> see that only +a has an influence, so latest proj version seems to use a
>> spherical version of eqc.
>
>
> If you look at
> <http://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_eqc.c>
>
> and the chapter 12 of Snyders manual, you will only find formulas for the
> sphere. So I guess there is no other way to calculate eqc.
>
> Maybe older versions calculated another radius for the sphere when an
> ellipsoid was given.

Stephen's "shift" was about 20km south, which correlates quite well if
you use semi-minor axis of WGS84 as radius of sphere while calculating
forward, and semi-major axis as radius of sphere while calculating the
inverse. At latitude of 55 degrees the difference is ca. 20 530 meters
(55 degrees -> 54.8156 degrees).

There are several different radii of the Earth, and some of them could
arguably be used in this context in place of semi-major axis.

All radii ar not created equal: http://en.wikipedia.org/wiki/Earth_radius

You have to bear in mind that this projection is intended for small
scale mapping, for example mapping the whole world. In that scale 20
km is nothing. If you need better representation of the Earth you have
to use a projection which takes ellipsoidal properties into account.
Of course the beef in this thread is not about choosing a projection,
but the change/difference in formulae used, which can create problems
as Stephen pointed out.

Cheers,

  Antti
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: [gdal-dev] gdalwarp EPSG:32662 problem

Stephen Woodbridge
On 2/10/2014 3:41 AM, Antti Castrén wrote:

> 2014-02-09 8:40 GMT+02:00 Andre Joost <[hidden email]>:
>> Am 09.02.2014 00:42, schrieb Even Rouault:
>>>
>>> But Antti guess seems right. Instead of +ellps=WGS84 (or +datum=WGS84), if
>>> you
>>> play with the a (semi-major axis) and b (semi-minor axis) parameters, you
>>> can
>>> see that only +a has an influence, so latest proj version seems to use a
>>> spherical version of eqc.
>>
>>
>> If you look at
>> <http://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_eqc.c>
>>
>> and the chapter 12 of Snyders manual, you will only find formulas for the
>> sphere. So I guess there is no other way to calculate eqc.
>>
>> Maybe older versions calculated another radius for the sphere when an
>> ellipsoid was given.
>
> Stephen's "shift" was about 20km south, which correlates quite well if
> you use semi-minor axis of WGS84 as radius of sphere while calculating
> forward, and semi-major axis as radius of sphere while calculating the
> inverse. At latitude of 55 degrees the difference is ca. 20 530 meters
> (55 degrees -> 54.8156 degrees).
>
> There are several different radii of the Earth, and some of them could
> arguably be used in this context in place of semi-major axis.
>
> All radii ar not created equal: http://en.wikipedia.org/wiki/Earth_radius
>
> You have to bear in mind that this projection is intended for small
> scale mapping, for example mapping the whole world. In that scale 20
> km is nothing. If you need better representation of the Earth you have
> to use a projection which takes ellipsoidal properties into account.
> Of course the beef in this thread is not about choosing a projection,
> but the change/difference in formulae used, which can create problems
> as Stephen pointed out.

Right, The real issue here what changed between Proj4 4.5 and Proj 4.8
that has caused this shift and how to best correct for it.

I hacked a solution into my problem by changing the reference bbox for
the image to shove it north by 20,000 meters. As you point out maybe I
should have pushed it 20,530 meters.

Anyway, thanks to everyone for taking the time to help me with this. It
is greatly appreciated.

-Steve
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj