[gdal-dev] problems using epsg:3785
Upendra Dadi
udadi at gmu.edu
Mon Nov 23 12:53:00 EST 2009
Peter,
One reason could be that your original file does not have EPSG:32622, but something else. Or it is defined incorrectly in gdal, though I doubt it. A quick comparison of EPSG:32622 in proj.csv and gdalinfo for the original file might give you some hint. Try warping original file to EPSG:32622(without using -s_srs) and and then do gdalinfo on the resulting file. Compare the resulting file with original, they should have same bbox if the original file has EPSG:32622. This is my one cent opinion.
Upendra
----- Original Message -----
From: Peter Freimuth <pfreimuth at arcor.de>
Date: Monday, November 23, 2009 11:30 am
Subject: [gdal-dev] problems using epsg:3785
> 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 --------------
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
More information about the gdal-dev
mailing list