[gdal-dev] Re: gdal to define projection, EPSG:3412

John Callahan john.callahan at UDel.Edu
Thu Nov 12 11:02:31 EST 2009 

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