[gdal-dev] Re: gdal to define projection, EPSG:3412
Hermann Peifer
peifer at gmx.eu
Thu Nov 12 12:43:23 EST 2009
John,
I guess you are right: "+lat_0=-90" shouldn't result into
"PARAMETER["latitude_of_origin",-70]". There seems to be some parameter
confusion.
It also looks like
http://spatialreference.org/ref/epsg/3412/prettywkt/ would be wrong in
this case.
EPSG guidance notes 7.2 describe polar stereographic projection with
variants A, B and C... I am not the ultimate expert here. It might be
worth to ask for clarification in the proj mailing list.
Hermann
John Callahan wrote:
> Thanks Frank/Hermann for your responses. However, something odd seems
> to be going on here. I feel like I'm missing something on how GDAL
> handles coordinate system information.
>
> Using the same input data set, and converting to different formats, the
> description of the output coordinate systems are different. Is this
> expected? Here are the gdal_translate commands and their respective
> output files gdalinfo descriptions below.
>
> [1] gdal_translate -a_srs "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0
> +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs" -of
> GTiff nt_200712_f13_v01_s.bsq nsidc_test.tif
>
> [2] gdal_translate -a_srs "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0
> +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs" -of HFA
> nt_200712_f13_v01_s.bsq nsidc_test.img
>
> [3] gdal_translate -a_srs "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0
> +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs" -of
> HDF4Image nt_200712_f13_v01_s.bsq nsidc_test.hdf
>
> [4] gdal_translate -a_srs EPSG:3412 -of GTiff nt_200712_f13_v01_s.bsq
> nsidc_test_epsg.tif
>
> [5] gdal_translate -a_srs EPSG:3412 -of HFA nt_200712_f13_v01_s.bsq
> nsidc_test_epsg.img
>
>
> Fortunately, I can display these images fine. They overlay nicely in
> QGIS with some vector data I have. However, the problem occurs when
> projecting to EPSG:4326 (or similar) coordinate system; the output files
> are incorrect/missing data. (It displays data from roughly lat -39.25
> only down to roughly -54.75 when data exists all the way down to -90.)
>
> One possible cause could be that gdalinfo reports (for all of these
> datasets) latitude_of_origin as -70. From
> http://nsidc.org/data/atlas/epsg_3412.html, the latitude of origin
> should be -90 and the latitude of standard parallel should be -70. Or
> the problem could be somewhere in gdalwarp, or specifically with the
> Hughes 1980 ellipsoid. I have no idea. Thanks.
>
> - John
>
>
>
> [1] Driver: GTiff/GeoTIFF
> Files: nsidc_test.tif
> Size is 316, 332
> Coordinate System is:
> PROJCS["unnamed",
> GEOGCS["unnamed ellipse",
> DATUM["unknown",
> SPHEROID["unnamed",6378273,298.279411123061]],
> PRIMEM["Greenwich",0],
> UNIT["degree",0.0174532925199433]],
> PROJECTION["Polar_Stereographic"],
> PARAMETER["latitude_of_origin",-70],
> PARAMETER["central_meridian",0],
> PARAMETER["scale_factor",1],
> PARAMETER["false_easting",0],
> PARAMETER["false_northing",0],
> UNIT["metre",1,
> AUTHORITY["EPSG","9001"]]]
> Origin = (-3950000.000000000000000,4350000.000000000000000)
> Pixel Size = (25000.000000000000000,-25000.000000000000000)
> Metadata:
> AREA_OR_POINT=Area
> Image Structure Metadata:
> INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left (-3950000.000, 4350000.000) ( 42d14'27.21"W, 39d13'51.20"S)
> Lower Left (-3950000.000,-3950000.000) (135d 0'0.00"W, 41d26'49.04"S)
> Upper Right ( 3950000.000, 4350000.000) ( 42d14'27.21"E, 39d13'51.20"S)
> Lower Right ( 3950000.000,-3950000.000) (135d 0'0.00"E, 41d26'49.04"S)
> Center ( 0.000, 200000.000) ( 0d 0'0.01"E, 88d 9'14.17"S)
> Band 1 Block=316x25 Type=Byte, ColorInterp=Gray
>
>
> [2] Driver: HFA/Erdas Imagine Images (.img)
> Files: nsidc_test.img
> Size is 316, 332
> Coordinate System is:
> PROJCS["Stereographic_South_Pole",
> GEOGCS["GCS_unnamed ellipse",
> DATUM["unknown",
> SPHEROID["Unknown",6378273,298.279411123061]],
> PRIMEM["Greenwich",0],
> UNIT["Degree",0.017453292519943295]],
> PROJECTION["Polar_Stereographic"],
> PARAMETER["latitude_of_origin",-70],
> PARAMETER["central_meridian",0],
> PARAMETER["false_easting",0],
> PARAMETER["false_northing",0],
> UNIT["Meter",1]]
> Origin = (-3950000.000000000000000,4350000.000000000000000)
> Pixel Size = (25000.000000000000000,-25000.000000000000000)
> Corner Coordinates:
> Upper Left (-3950000.000, 4350000.000) ( 42d14'27.21"W, 39d13'51.20"S)
> Lower Left (-3950000.000,-3950000.000) (135d 0'0.00"W, 41d26'49.04"S)
> Upper Right ( 3950000.000, 4350000.000) ( 42d14'27.21"E, 39d13'51.20"S)
> Lower Right ( 3950000.000,-3950000.000) (135d 0'0.00"E, 41d26'49.04"S)
> Center ( 0.000, 200000.000) ( 0d 0'0.01"E, 88d 9'14.17"S)
> Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined
> Metadata:
> LAYER_TYPE=athematic
>
>
> [3] Driver: HDF4Image/HDF4 Dataset
> Files: nsidc_test.hdf
> Size is 316, 332
> Coordinate System is:
> PROJCS["unnamed",
> GEOGCS["unnamed ellipse",
> DATUM["unknown",
> SPHEROID["unnamed",6378273,298.279411123061]],
> PRIMEM["Greenwich",0],
> UNIT["degree",0.0174532925199433]],
> PROJECTION["Polar_Stereographic"],
> PARAMETER["latitude_of_origin",-70],
> PARAMETER["central_meridian",0],
> PARAMETER["scale_factor",1],
> PARAMETER["false_easting",0],
> PARAMETER["false_northing",0],
> UNIT["Meter",1]]
> Origin = (-3950000.000000000000000,4350000.000000000000000)
> Pixel Size = (25000.000000000000000,-25000.000000000000000)
> Metadata:
> Signature=Created with GDAL (http://www.remotesensing.org/gdal/)
> TransformationMatrix=-3950000.000000, 25000.000000, 0.000000,
> 4350000.000000, 0.000000, -25000.000000
> Projection=PROJCS["unnamed",GEOGCS["unnamed
> ellipse",DATUM["unknown",SPHEROID["unnamed",6378273,298.279411123061]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
>
> Corner Coordinates:
> Upper Left (-3950000.000, 4350000.000) ( 42d14'27.21"W, 39d13'51.20"S)
> Lower Left (-3950000.000,-3950000.000) (135d 0'0.00"W, 41d26'49.04"S)
> Upper Right ( 3950000.000, 4350000.000) ( 42d14'27.21"E, 39d13'51.20"S)
> Lower Right ( 3950000.000,-3950000.000) (135d 0'0.00"E, 41d26'49.04"S)
> Center ( 0.000, 200000.000) ( 0d 0'0.01"E, 88d 9'14.17"S)
> Band 1 Block=316x316 Type=Byte, ColorInterp=Gray
>
>
> [4] Driver: GTiff/GeoTIFF
> Files: nsidc_test_epsg.tif
> Size is 316, 332
> Coordinate System is:
> PROJCS["NSIDC Sea Ice Polar Stereographic South",
> GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",
> DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",
> SPHEROID["Hughes 1980",6378273,298.279411123061,
> AUTHORITY["EPSG","7058"]],
> AUTHORITY["EPSG","6054"]],
> PRIMEM["Greenwich",0],
> UNIT["degree",0.0174532925199433],
> AUTHORITY["EPSG","4054"]],
> UNIT["metre",1,
> AUTHORITY["EPSG","9001"]],
> AUTHORITY["EPSG","3412"]]
> Origin = (-3950000.000000000000000,4350000.000000000000000)
> Pixel Size = (25000.000000000000000,-25000.000000000000000)
> Metadata:
> AREA_OR_POINT=Area
> Image Structure Metadata:
> INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left (-3950000.000, 4350000.000)
> Lower Left (-3950000.000,-3950000.000)
> Upper Right ( 3950000.000, 4350000.000)
> Lower Right ( 3950000.000,-3950000.000)
> Center ( 0.000, 200000.000)
> Band 1 Block=316x25 Type=Byte, ColorInterp=Gray
>
>
> [5] Driver: HFA/Erdas Imagine Images (.img)
> Files: nsidc_test_epsg.img
> Size is 316, 332
> Coordinate System is:
> PROJCS["NSIDC_Sea_Ice_Polar_Stereographic_South",
> GEOGCS["GCS_Unspecified datum based upon the Hughes 1980 ellipsoid",
> DATUM["Hungarian_Datum_1909",
> SPHEROID["Hughes_1980",6378273,298.279411123061]],
> PRIMEM["Greenwich",0],
> UNIT["Degree",0.017453292519943295]],
> PROJECTION["Polar_Stereographic"],
> PARAMETER["latitude_of_origin",-70],
> PARAMETER["central_meridian",0],
> PARAMETER["false_easting",0],
> PARAMETER["false_northing",0],
> UNIT["Meter",1]]
> Origin = (-3950000.000000000000000,4350000.000000000000000)
> Pixel Size = (25000.000000000000000,-25000.000000000000000)
> Corner Coordinates:
> Upper Left (-3950000.000, 4350000.000) ( 42d14'27.21"W, 39d13'51.20"S)
> Lower Left (-3950000.000,-3950000.000) (135d 0'0.00"W, 41d26'49.04"S)
> Upper Right ( 3950000.000, 4350000.000) ( 42d14'27.21"E, 39d13'51.20"S)
> Lower Right ( 3950000.000,-3950000.000) (135d 0'0.00"E, 41d26'49.04"S)
> Center ( 0.000, 200000.000) ( 0d 0'0.01"E, 88d 9'14.17"S)
> Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined
> Metadata:
> LAYER_TYPE=athematic
>
>
>
>
> Frank Warmerdam wrote:
>> Hermann Peifer wrote:
>>> John Callahan wrote:
>>>> Thanks for responding Hermann. Yes, I tried using -s_srs EPSG:3412
>>>> as well as the full PROJ4 string with the same result. When
>>>> specifying -s_srs, does gdalwarp ignore (override) the project
>>>> information in the hdr file?
>>>>
>>>
>>> Hmm. AFAIK: gdal_translate -a_srs does override the file's projection
>>> information. I have a random.tif in LAEA projection (EPSG:3035), so I
>>> tried the following with gdal 1.6.2:
>>>
>>> gdal_translate random.tif out1.tif -a_srs EPSG:3412 # this results
>>> in [1]
>>>
>>> gdal_translate random.tif out2.tif -a_srs "+proj=stere +lat_0=-90
>>> +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449
>>> +units=m +no_defs" # this gives [2]
>>>
>>> To me as an innocent GDAL user, this behaviour looks similar to
>>> http://trac.osgeo.org/gdal/ticket/3016. I might be wrong, though, in
>>> particular as you wrote that using the PROJ.4 string didn't help (and
>>> your source file isn't a GeoTIFF in the first place).
>>
>> Hermann,
>>
>> I have confirmed this problem and filed it as:
>>
>> http://trac.osgeo.org/gdal/ticket/3220
>>
>> It is presumably a problem with translation of a specific projection
>> method,
>> possibly in libgeotiff itself.
>>
>> Best regards,
More information about the gdal-dev
mailing list