[postgis-users] [Fossgis-talk] problems using epsg:3785
dassouki
dassouki at gmail.com
Mon Nov 23 08:27:50 PST 2009
Hopefully my bad day is over now that i spilled yogurt on my lap lool. i
laughed. So how's the pomegranate ? how's your day going so far?
Jan Hartmann wrote:
> 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
> ------------------------------------------------------------------------
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
More information about the postgis-users
mailing list