[gdal-dev] problems using epsg:3785
Peter Freimuth
pfreimuth at arcor.de
Mon Nov 23 11:30:02 EST 2009
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 --------------
Driver: GTiff/GeoTIFF
Files: 2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB_pix-3785.tif
Size is 5051, 5081
Coordinate System is:
PROJCS["unnamed",
GEOGCS[,
DATUM["unknown",
SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3785"]]
Origin = (-5470299.855895809800000,-2551883.714705614400000)
Pixel Size = (5.421926883630149,-5.421926883630149)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-5470299.856,-2551883.715)
Lower Left (-5470299.856,-2579432.525)
Upper Right (-5442913.703,-2551883.715)
Lower Right (-5442913.703,-2579432.525)
Center (-5456606.780,-2565658.120)
Band 1 Block=5051x1 Type=Byte, ColorInterp=Red
Band 2 Block=5051x1 Type=Byte, ColorInterp=Green
Band 3 Block=5051x1 Type=Byte, ColorInterp=Blue
-------------- next part --------------
Driver: GTiff/GeoTIFF
Files: 2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB_32622-pix-3785.tif
Size is 5066, 5066
Coordinate System is:
PROJCS["unnamed",
GEOGCS[,
DATUM["unknown",
SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3785"]]
Origin = (-5470299.855895809800000,-2535649.341730865200000)
Pixel Size = (5.406746173490638,-5.406746173490638)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-5470299.856,-2535649.342)
Lower Left (-5470299.856,-2563039.918)
Upper Right (-5442909.280,-2535649.342)
Lower Right (-5442909.280,-2563039.918)
Center (-5456604.568,-2549344.630)
Band 1 Block=5066x1 Type=Byte, ColorInterp=Red
Band 2 Block=5066x1 Type=Byte, ColorInterp=Green
Band 3 Block=5066x1 Type=Byte, ColorInterp=Blue
-------------- next part --------------
Driver: PCIDSK/PCIDSK Database File
Files: 2009-03-02T140133_RE3_3A-NAC_642351_34785_enhancedRGB.pix
Size is 5000, 5000
Coordinate System is:
PROJCS["UTM Zone 22, Northern Hemisphere",
GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",
DATUM["Not specified (based on WGS 84 spheroid)",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-51],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (691500.000000000000000,-2471500.000000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Corner Coordinates:
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)
Center ( 704000.000,-2484000.000) ( 49d 1'3.43"W, 22d27'0.49"S)
Band 1 Block=5000x1 Type=Byte, ColorInterp=Undefined
Description = Contents Not Specified
...
Band 2 Block=5000x1 Type=Byte, ColorInterp=Undefined
Description = Contents Not Specified
...
Band 3 Block=5000x1 Type=Byte, ColorInterp=Undefined
Description = Contents Not Specified
...
More information about the gdal-dev
mailing list