[postgis-users] [Fossgis-talk] problems using epsg:3785
Jan Hartmann
j.l.h.hartmann at uva.nl
Mon Nov 23 08:16:52 PST 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20091123/97ccfd69/attachment.html>
More information about the postgis-users
mailing list