[postgis-users] problems using epsg:3785
Peter Freimuth
pfreimuth at arcor.de
Mon Nov 23 07:16:20 PST 2009
Hi all,
i am struggling around with a projection problem. I want to harmonize
several hundreds of satellite images from different UTM Zones to the
Google Projection Spherical Mercator. I am using a python script which
handles the job for me using "gdalwarp" from the 1.6.1 and 1.7.0dev
release (same result).
Find attached the gdalinfo output of the original file
(gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB.txt):
Upper Left ( 691500.000,-2471500.000) ( 49d 8'25.94"W, 22d20'19.37"S)
Lower Left ( 691500.000,-2496500.000) ( 49d 8'15.09"W, 22d33'52.01"S)
Upper Right ( 716500.000,-2471500.000) ( 48d53'52.47"W, 22d20'8.69"S)
Lower Right ( 716500.000,-2496500.000) ( 48d53'40.21"W, 22d33'41.21"S)
echo "warp giving a source epsg code"
gdalwarp -s_srs "EPSG:32622" -t_srs "EPSG:3785"
2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB.pix
c:\tmp\2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB_32622-pix-3785.tif
this results in
(gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_32622-pix-3785.txt)
Upper Left (-5470299.856,-2535649.342)
Lower Left (-5470299.856,-2563039.918)
Upper Right (-5442909.280,-2535649.342)
Lower Right (-5442909.280,-2563039.918)
echo "warp without giving a source epsg code"
gdalwarp -t_srs "EPSG:3785"
2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB.pix
c:\tmp\2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB_pix-3785.tif
This results in
(gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_pix-3785.txt)
Upper Left (-5470299.856,-2551883.715)
Lower Left (-5470299.856,-2579432.525)
Upper Right (-5442913.703,-2551883.715)
Lower Right (-5442913.703,-2579432.525)
Somehow i get big shift in the Y values of the images.
As a kind of reference i have a polygon layer in epsg:4326 stored in a
postgis 1.4.0 database table which contains the exact definitions of the
tiling grid which is used to create the satallite images.
The same shift can be reproduce in postgis when using a somehow
different epsg:3785 definition:
As you can see, there is a significant difference in the Y value.
select
astext(envelope(transform(GeomFromEWKT('SRID=4326;POLYGON((-49.1405410766602
-22.5644474029541,-49.1405410766602 -22.3357467651367,-48.8945007324219
-22.3357467651367,-48.8945007324219 -22.5644474029541,-49.1405410766602
-22.5644474029541))'::TEXT),3785)))
=>>"POLYGON((-5470300 -2579430.25,-5470300 -2551883.5,-5442911
-2551883.5,-5442911 -2579430.25,-5470300 -2579430.25))"
select
astext(envelope(transform(GeomFromEWKT('SRID=4326;POLYGON((-49.1405410766602
-22.5644474029541,-49.1405410766602 -22.3357467651367,-48.8945007324219
-22.3357467651367,-48.8945007324219 -22.5644474029541,-49.1405410766602
-22.5644474029541))'::TEXT),903785)))
=>>"POLYGON((-5470300 -2563038,-5470300 -2535649.25,-5442911
-2535649.25,-5442911 -2563038,-5470300 -2563038))"
Here are the two definitons as they are in my PostGIS spatial_ref_sys table.
"3785" as defined in postgis 1.4.0 and proj-4.7.0 epsg file
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+units=m +k=1.0 +nadgrids=@null +no_defs
PROJCS["Popular Visualisation CRS / Mercator (deprecated)",
GEOGCS["Popular Visualisation CRS",
DATUM["Popular_Visualisation_Datum",
SPHEROID["Popular Visualisation Sphere",6378137,0,
AUTHORITY["EPSG","7059"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6055"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4055"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","3785"],
AXIS["X",EAST],
AXIS["Y",NORTH]]"
I added a second definition which i took from spatialreference.org
(http://spatialreference.org/ref/epsg/3785) with the srid 903785
"903785"
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs
PROJCS["Popular Visualisation CRS / Mercator",
GEOGCS["Popular Visualisation CRS",
DATUM["Popular_Visualisation_Datum",
SPHEROID["Popular Visualisation Sphere",6378137,0,
AUTHORITY["EPSG","7059"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6055"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4055"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","3785"],
AXIS["X",EAST],
AXIS["Y",NORTH]]
Now the big questions:
Which of the above given definitions is the correct one?? So which of
the polygons is the correct result when i want epsg:3785?
How can it be explained that gdalwarp seems to use two different
projections when warping the same file to the same output projection?
Any ideas or comments??
Thanks for any kind of help,
Peter
--
Peter Freimuth
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_pix-3785.txt
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20091123/c79e6331/attachment.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_32622-pix-3785.txt
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20091123/c79e6331/attachment-0001.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: gdalinfo_2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB.txt
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20091123/c79e6331/attachment-0002.txt>
More information about the postgis-users
mailing list