[Xastir] DRG Problem

[Tom Russo]

On Fri, Mar 11, 2005 at 01:58:18PM -0700, we recorded a bogon-computron collision of the <russo at bogodyn.org> flavor, containing:
> On Fri, Mar 11, 2005 at 01:03:37PM -0700, we recorded a bogon-computron collision of the <russo at bogodyn.org> flavor, containing:
> > On Fri, Mar 11, 2005 at 02:39:21PM -0500, we recorded a bogon-computron collision of the <mckeehan at mckeehan.homeip.net> flavor, containing:
> > > One more piece of information:
> > > 
> > > The Digital Raster Graphic (DRG) is a digital raster image of a scanned U.S.
> > > Geological Survey (USGS) topographic map georeferenced to the State Plane
> > > Coordinate System 1983 grid.
> > 
> > AH.  If that's the Tennessee state plane system of 1983, it's EPSG number
> > 32136.  
> > 
> > But using the .tfw file leads to some odd information.
> 
> I've been beating on it for a bit, and am no closer.  

Actually, I think I see.  I almost had it.  The big problem is that TVA's
version of SPCS 83 uses US survey feet.  I tried to play with proj.4 parameters
to get a valid cs2cs parameter set that used US survey feet, but 
For whatever reason, my cs2cs command with US-Feet isn't doing the right thing.
I must be missing some subtlety of how it uses its "+units" argument.
No matter what I tried, I could not get a command that would convert coordinates
in feet given by TVA into the corresponding lat/lon points.

But if I do all lat/lon<->state plane coordinates in meters and do 
manual conversions between feet and meters, things seem alright.
  (i.e. if I do cs2cs -v +init="epsg:32136" +to +proj=latlong +datum=NAD83
   and give it the coordinates I get from the tfw files after converting
   feet to meters,  it gives sane values of latitude and longitude)

So here's what I did:

  awk '{print $1*.3048}' < C35084A1.tfw > C35084A1_meter.tfw
This converts all the parameters in the tfw file to meters instead of feet.
The TFW file defines the linear affine transformation between pixel coordinates
and the projected coordinate system.  So once this new tfw file is attached to
the tif file, one can easily determine the SPCS coordinates IN METERS of
each pixel.

  geotifcp -e C35084A1_meter.tfw C35084A1.tif junk.tif
This tells geotifcp to install the positioning info if the .tfw file into
the TIF image.  This step gives me a partial GeoTIFF with coordinates in 
meters, as you can now see by doing a "listgeo" on the junk.tif file.  Since 
the projection isn't yet stored in GeoTIFF tags, all listgeo can show you is 
the SPCS coordinates of the corners and center.

  gdal_translate -a_srs EPSG:32136 junk.tif C35084A1_meter.tif
This puts in all the metadata necessary to show that C35084A1_meter.tif is
georeferenced to the state plane coordinate system for Tennessee IN METERS
(which is what EPSG defines it as, even though TVA uses survey feet).  This
only works because the image is already in those coordinates --- we're just
adding the GeoTIFF tags necessary to let programs figure out what to do
with them.  Listgeo will now show that this new file has tons of new tags 
in it.  Since it knows the projection now, it can also convert the corner 
coordinates to lat/lon, and prints them out alongside the SPCS coordinates.  
They look right to me.

Once you're there, you can use gdalwarp to change its coordinate system.
Ideally, one should be able to do:

 gdalwarp -t_srs EPSG:4326 C35084A1_meter.tif C35084A1_latlon.tif

to get the map warped into lat/lon, which will be the best you can get out of 
xastir.  I've done this for all sorts of Geotiff images, BUT for some odd 
reason this isn't working out for me on this set of data --- it gives 
bogus coordinates for the corner points.  Your milage may vary --- I am 
running gdal out of CVS, and perhaps there's something broken.  I tried to
run gdalwarp on other files I've converted before, and it worked.  I don't
know how to get around this.

HOWEVER, I *can* do:
  gdalwarp -t_srs EPSG:26916 C35084A1_meter.tif C35084A1_utm.tif
to get a UTM version of the map whose coordinates ARE correct.  If I attempt
to warp *this* to lat/lon it also fails (leading me to question whether my
gdal is broken).  But xastir can (approximately) handle utm geotiffs so
this could be OK for you.  Unfortunately, xastir's approximate handling of 
UTM is optimized (quite well) for 7.5' (1:24000 scale) quads, and on something 
as big as the 1:250000 there's noticable distortion and incorrect positioning 
(but not as bad as the .geo stuff would be).

I still think you're better off hunting for properly georeferenced GeoTIFFS
rather than trying to hack on the TVA data like this.  But if TVA is the
best you can do, you can follow this procedure to get them into something
xastir can handle.  Your best bet would be to use the 7.5' quads you can get
from TVA, because the distortion will be minimum when displayed in xastir
as UTM maps.  If you can get the "gdalwarp -t_srs EPSG:4326" version to work,
then you can use any of the maps because that'll convert the maps directly
to the coordinate system xastir uses and xastir's approximate transformations
won't be an issue.

-- 
Tom Russo    KM5VY     SAR502  DM64ux         http://www.swcp.com/~russo/
Tijeras, NM  QRPL#1592 K2#398  SOC#236 AHTB#1 
 "When life gives you lemons, find someone with a paper cut."