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

Thu, 12 Nov 2009 13:28:26 -0800

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
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev






















					Reply via email to