[Gdal-dev] PROJ.4 warping string ignored
Seth Price
seth at pricepages.org
Thu Nov 1 14:40:39 EDT 2007
Thanks for your help. I think I might have a bug for you. Using the
file I've produced with your hack ('+proj=merc +a=6378137.0
+b=6378137.0 +nadgrids=@null +wktext +units=m') I was able to produce
the expected output file.
However, when I attempt to break that output into tiles (using
'gdalwarp' and the same -t_srs string) the output image is shifted (as
if we are back in the ellipsoid I think). However, the projection info
from the file remains identical. I found that removing the 'nadgrids'
hack produces the expected behavior.
I thought I'd report this bug in case it is, in fact, a bug. Here is
the output from the command line that demonstrates what I mention
above. 1_03868.no_srs.tif and 1_03868.srs.tif are images of different
areas.
$ ./gdalinfo N010W074_x1.merc.tif
Driver: GTiff/GeoTIFFSize is 11651, 11797
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6378137,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-8237642.318702000193298,1118889.974857999943197)
Pixel Size = (19.109001938546044,-19.108775438247005)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=219
TIFFTAG_YRESOLUTION=222
TIFFTAG_RESOLUTIONUNIT=3 (pixels/cm)
Corner Coordinates:
Upper Left (-8237642.319, 1118889.975) ( 74d 0'0.00"W, 10d 0'0.00"N)
Lower Left (-8237642.319, 893463.751) ( 74d 0'0.00"W, 8d 0'0.00"N)
Upper Right (-8015003.337, 1118889.975) ( 72d 0'0.00"W, 10d 0'0.00"N)
Lower Right (-8015003.337, 893463.751) ( 72d 0'0.00"W, 8d 0'0.00"N)
Center (-8126322.828, 1006176.863) ( 73d 0'0.00"W, 9d 0'4.98"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Overviews: 3884x3933, 1942x1967, 971x984, 486x492, 243x246, 122x123
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Overviews: 3884x3933, 1942x1967, 971x984, 486x492, 243x246, 122x123
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
Overviews: 3884x3933, 1942x1967, 971x984, 486x492, 243x246, 122x123
$ ./gdalwarp -of GTiff -ts 256 256 -te -8233185.1903991
1110477.1468928 -8228293.220589 1115369.1167029 N010W074_x1.merc.tif
1_03868.no_srs.tifCreating output file that is 256P x 256L.
Processing input file N010W074_x1.merc.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
$ ./gdalinfo 1_03868.no_srs.tif
Driver: GTiff/GeoTIFFSize is 256, 256
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6378137,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-8233185.190399100072682,1115369.116702900035307)
Pixel Size = (19.109257070704189,-19.109257070703279)
Metadata:
AREA_OR_POINT=Area
Corner Coordinates:
Upper Left (-8233185.190, 1115369.117) ( 73d57'35.86"W, 9d58'7.86"N)
Lower Left (-8233185.190, 1110477.147) ( 73d57'35.86"W, 9d55'32.04"N)
Upper Right (-8228293.221, 1115369.117) ( 73d54'57.66"W, 9d58'7.86"N)
Lower Right (-8228293.221, 1110477.147) ( 73d54'57.66"W, 9d55'32.04"N)
Center (-8230739.205, 1112923.132) ( 73d56'16.76"W, 9d56'49.95"N)
Band 1 Block=256x32 Type=Byte, ColorInterp=Red
Band 2 Block=256x32 Type=Byte, ColorInterp=Green
Band 3 Block=256x32 Type=Byte, ColorInterp=Blue
$ ./gdalwarp -of GTiff -ts 256 256 -t_srs '+proj=merc +a=6378137.0
+b=6378137.0 +nadgrids=@null +wktext +units=m' -te -8233185.1903991
1110477.1468928 -8228293.220589 1115369.1167029 N010W074_x1.merc.tif
1_03868.srs.tif
Creating output file that is 256P x 256L.
Processing input file N010W074_x1.merc.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
$ ./gdalinfo 1_03868.srs.tif
Driver: GTiff/GeoTIFFSize is 256, 256
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6378137,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-8233185.190399100072682,1115369.116702900035307)
Pixel Size = (19.109257070704189,-19.109257070703279)
Metadata:
AREA_OR_POINT=Area
Corner Coordinates:
Upper Left (-8233185.190, 1115369.117) ( 73d57'35.86"W, 9d58'7.86"N)
Lower Left (-8233185.190, 1110477.147) ( 73d57'35.86"W, 9d55'32.04"N)
Upper Right (-8228293.221, 1115369.117) ( 73d54'57.66"W, 9d58'7.86"N)
Lower Right (-8228293.221, 1110477.147) ( 73d54'57.66"W, 9d55'32.04"N)
Center (-8230739.205, 1112923.132) ( 73d56'16.76"W, 9d56'49.95"N)
Band 1 Block=256x32 Type=Byte, ColorInterp=Red
Band 2 Block=256x32 Type=Byte, ColorInterp=Green
Band 3 Block=256x32 Type=Byte, ColorInterp=Blue
On Oct 24, 2007, at 5:05 PM, Mateusz Loskot wrote:
> Seth,
>
>
> I found that your command is wrong.
> The OGR is correct, but you need to use +nadgrids=@null trick to
> overwrite ellipsoid issue.
>
>
> It is well explained in PROJ.4 FAQ:
>
>
> http://proj.maptools.org/faq.html (Changing Ellipsoid...)
>
>
>
>
>
>
>
>
> Seth Price wrote:
>> Sorry if I'm being impatient. Any news on this bug?
>> Thanks,
>> Seth
>>
>>
>> On Oct 23, 2007, at 9:18 AM, Mateusz Loskot wrote:
>>
>>> Seth Price wrote:
>>>> Thanks for the response. Although I don't claim to know what I'm
>>>> doing,
>>>> I'm trying to define a sphere (not the normal ellipsoid).
>>>
>>> Seth,
>>>
>>> Ah, right. I've missed the b semi-axis parameter.
>>>
>>>> I'm copying
>>>> the instructions from this page, the PROJ.4 parameters are at the
>>>> bottom
>>>> in the 'Update 4' text:
>>>> http://cfis.savagexi.com/articles/2006/05/03/google-maps-deconstructed
>>>>
>>>> Also, it doesn't make sense that I'm defining the same parameters,
>>>> because if I directly use the proj utility, the output is
>>>> different:
>>>>
>>>> proj +proj=merc +datum=WGS84 +units=m -r << EOF
>>>> ? 46 -90
>>>> ? EOF
>>>> -10018754.17 5749599.55
>>>>
>>>> proj +proj=merc +latts=0 +lon0=0 +k=1.0 +x0=0 +y0=0 +a=6378137.0
>>>> +b=6378137.0 +units=m -r << EOF
>>>> ? 46 -90
>>>> ? EOF
>>>> -10018754.17 5780349.22
>>>>
>>>> I think that gdalwarp is confusing the two commands because the
>>>> image is
>>>> shifted the same amount as the difference of these two coordinates.
>>>
>>> Yes, your assumptions seem to be correct.
>>>
>>> I've checked the code, and construction of spatial reference system
>>> object (OGRSpatialReference class) ignores your a and b parameters,
>>> recognizes proj=merc name and uses default ellipsoid.
>>> IOW, according to my understanding, handling custom PROJ.4
>>> definition
>>> that uses standard projection name is not supported.
>>>
>>> Let's wait for Frank's comment.
>>>
>>> Cheers
>>> --
>>> Mateusz Loskot
>>> http://mateusz.loskot.net
>>
>>
>
>
>
>
> --
> Mateusz Loskot
> http://mateusz.loskot.net
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/gdal-dev
More information about the gdal-dev
mailing list