[Fossgis-talk] problems using epsg:3785

Jan Hartmann j.l.h.hartmann at uva.nl
Mo Nov 23 17:16:52 CET 2009


Hi Peter,

The two projection definitions you give for Google Mercator differ only 
in the datum shift:

+nadgrids=@null (epsg and PostGIS)
+towgs84=0,0,0,0,0,0,0 (spatialreference.org)

Number one works OK for me, number two is  wrong. You can find out about 
nadgrids at:

http://trac.osgeo.org/proj/wiki/FAQ, section "Changing Ellipsoid / Why 
can't I convert from WGS84 to Virtual Globe Mercator?"

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?

Jan

Peter Freimuth wrote:
> 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
>
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> Fossgis-talk-liste mailing list
> Fossgis-talk-liste at fossgis.de
> https://lists.fossgis.de/mailman/listinfo/fossgis-talk-liste
-------------- nächster Teil --------------
Ein Dateianhang mit HTML-Daten wurde abgetrennt...
URL: <https://lists.fossgis.de/pipermail/fossgis-talk-liste/attachments/20091123/97ccfd69/attachment-0002.htm>