<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html;charset=ISO-8859-15"
http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
Hi Peter,<br>
<br>
The two projection definitions you give for Google Mercator differ only
in the datum shift:<br>
<br>
+nadgrids=@null (epsg and PostGIS)<br>
+towgs84=0,0,0,0,0,0,0 (spatialreference.org)<br>
<br>
Number one works OK for me, number two is wrong. You can find out
about nadgrids at:<br>
<br>
<a class="moz-txt-link-freetext" href="http://trac.osgeo.org/proj/wiki/FAQ">http://trac.osgeo.org/proj/wiki/FAQ</a>, section "Changing Ellipsoid / Why
can't I convert from WGS84 to Virtual Globe Mercator?"<br>
<br>
As to the two gdalwarp's: *both* projection definitions need a
+datum, +towgs84 or a +nadgrids parameter, else the
WGS84-datum-transformation won't work. In your case, the first warp
looks OK to me (both projections are OK), the second one gets its
source projection from the source file. I would guess that something is
wrong with the datum information there, and that the first warp is the
correct one. Don't you have any ground information?<br>
<br>
Jan<br>
<br>
Peter Freimuth wrote:
<blockquote cite="mid:4B0AA744.1070202@arcor.de" type="cite">
<pre wrap="">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
(<a class="moz-txt-link-freetext" href="http://spatialreference.org/ref/epsg/3785">http://spatialreference.org/ref/epsg/3785</a>) 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
</pre>
<pre wrap=""><pre wrap=""><pre wrap=""><pre wrap="">
<hr size="4" width="90%">
_______________________________________________
Fossgis-talk-liste mailing list
<a class="moz-txt-link-abbreviated" href="mailto:Fossgis-talk-liste@fossgis.de">Fossgis-talk-liste@fossgis.de</a>
<a class="moz-txt-link-freetext" href="https://lists.fossgis.de/mailman/listinfo/fossgis-talk-liste">https://lists.fossgis.de/mailman/listinfo/fossgis-talk-liste</a>
</pre></pre></pre></pre>
</blockquote>
</body>
</html>