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 <[email protected]> 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 +nadgri...@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 > > >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
