[gdal-dev] Re: gdal to define projection, EPSG:3412
John Callahan
john.callahan at UDel.Edu
Thu Nov 12 16:28:07 EST 2009
I just received word back from NSIDC. They mentioned this is a known
issue when projecting from polar to cylindrical coordinate systems.
http://trac.osgeo.org/gdal/ticket/2729
The work around is to specify the output extents in the command. For my
data, the proper command is:
gdalwarp -of GTiff -s_srs EPSG:3412 -t_srs EPSG:4326 -te -180.0 -90 180.0 -39.2308884002772 nt_200712_f13_v01_s.bsq out.tif
Using -s_srs EPSG:3412 or the full PROJ.4 string didn't matter. Both
worked great.
This still doesn't explain why "+lat_0=-90" results into
"PARAMETER["latitude_of_origin",-70]" in the description. Since
gdalwarp does indeed work (with specifying the output extents), maybe
the error is in how gdalinfo returns the srs info (differently in
different formats, as Frank explained) and not somewhere inside gdalwarp
projection algorithms.
- John
Hermann Peifer wrote:
> 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,
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
More information about the gdal-dev
mailing list