[Gdal-dev] Re: Another GeoTif question

[Frank Warmerdam]

Stephen H. Savage wrote:
> Hi Frank,
> 
>    I've got a couple more questions for you about the geotifcp and gdalwarp 
> programs.  The first one requires a bit of explanation.
> 
> 1)  Part of the header information from a geo tiff file projected in WGS 84 UTM 
> zone 36N reads as follows:
> 
> ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_36N
> PCSCitationGeoKey (Ascii,24): "UTM International WGS84"
> 
> pcs.csv in GDAL\data has  32636	WGS 84 / UTM zone 36N, which comes from EPSG.
> 
> MapObjects has:	moProjCS_WGS1984UTM_36N	32636	WGS 1984 UTM Zone 36N
> 
>    Each of these formats is just a bit different.  My question involves trying 
> to understand just how sensitive the .tiff files are to specific syntax in the 
> header.  I would like to be able to extract the information from the MapObjects 
> data, because it is easily bound to a combo box control using a strings 
> collection in VB, thus eliminating the need to find the records in the excel 
> spreadsheets.  Can I use the MapObjects string directly to build the .lgo file, 
> or do I have to somehow duplicate exactly what the geotiff header has in it for 
> any projection?  What would be really nice would be to be able to plug the 
> numeric code for the projection into the .lgo file, rather than have to format 
> the whole set of entries.  Is there a way to do that?

Stephen,

The strings associated with the ProjectedCSTypeGeoKey must match the internal
list in geotifcp exactly in order for them to work when reading a metadata
file.  The PCSCitationGeoKey contents are never interpreted by geotifcp, and
are freeform text.  In short, you are not going to be able to trivally use
the mapobjects text to generate the geotiff, file, but the "32636" is the
actual numerical value that gets assigned to the ProjectedCSTypeGeoKey
field in the GeoTIFF file.  If you need to operate through a geotiff metadata
file then you can lookup the number to get the corresponding symbol using
the epsg_pcs.inc file distributed in libgeotiff.

>    My second question has to do with the gdalwarp program.  I can easily 
> extract the EPSG numeric values for the projection I want to shift to, and send 
> it with the EPSG:32636 format.  Does gdalwarp perform datum transformations 
> automatically, or do I have to do that explicitly somehow?  (For example, 
> reprojecting an image in UTM NAD 27, zone 12 to utm NAD 83 zone 12?).  Does the 
> program handle a two step transformation, for example, from the Palestine Belt 
> to WGS 84 UTM zone 36N - where I have to go from the Pal Belt to WGS 84 decimal 
> degrees first, then to the UTM projection?

gdalwarp will apply datum shifts in situations where it can figure out which
one is appropriate.  Note that NAD27 and NAD83 do not apply outside the Americas
to the best of my knowledge.  If a coordinate system is determined to be NAD27
gdalwarp will attempt to use the conus or ntv1_can.dat files to do the
NAD27/NAD83 shift.  Other areas (like Alaska, Hawaii) will be difficult.

It will also support "two step transformations".  You can give two distinct
PCSes and internally it will go through various steps, including an intermediate
WGS84 decimal degrees step.

Note that there are many situations in which gdalwarp cannot figure out the
appropriate datum shift to use.  Generally speaking there are just a few cases.

1) NAD27 - uses the grid shift files to convert to NAD83.
2) NAD83 - treated as identical to WGS84.
3) WGS84 - the reference datum.
4) TOWGS84[3parms] - systems with a TOWGS84[] clause can be automatically
            converted to WGS84.
5) TOWGS84[7parms] - systems with a TOWGS84[] clause can be automatically
            converted to WGS84.

The code that creates a coordinate system from an EPSG code will *not* produce
a TOWGS84[] clause if there are more than one possible approximation available.
Thus, things that seem like they should be transformable will often not be
because the system doesn't know what approximation to use.

>    I've got my little program reading MrSid or .tif files, and getting the 
> headers.  Then I can crop either of these two formats and save to .tif.  It 
> will currently write revised geotiff headers based on UTM zones, but not 
> anything else.  I want to support a wider range of projected and geographic 
> coordinate systems, and use the gdalwarp program to reproject them from the 
> system read from their current header to whatever I select from a drop down 
> list.  After that, I plan to implement the gdal_merge.py routine.  
> 
>    If I can get it all to work, I'll probably bundle it with an expanded 
> version of my ReprojectME! program, which currently reprojects ArcView 
> shapefiles, ArcInfo coverages, ascii coordinate lists and individual points 
> between about 24 projections used in the Middle East.  The program can be 
> downloaded free from my web page, by pressing the "ReprojectME! Program" 
> picture on the opening screen and following the links from there.  It would be 
> really spiffy to have a little utility that will handle vector, raster and 
> ascii data.  Something to keep me busy for a bit...

Very cool!

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent