
Greetings.  I'm using grass5beta11 on my FreeBSD 4.3-stable machine.  I'm
trying to use USGS Digital Line Graphs for my area to play around with the
system and learn more about it.  I just subscribed to GRASSLIST today, but
I've perused the archives and don't see the questions I have below answered.

Because I'm most familiar with United States Geological Survey 7.5' quad maps,
I'm using the USGS "SDTS" format DEM and DLG files available for free download
(see below, in the gorey details section).  As I play with Grass, I'm able to
compare what I see on my screen with the USGS quad map of the same area I have
on paper in front of me.  It's been fun, but I've hit a wall in importing some
data.

My message has three parts.  A simple question, a slightly detailed
restatement of the same question, and a long-winded, overly detailed
restatement of how I came to conclude that I have a problem.

Simple statement of the question: How can I import two DLGs that were created
using different geodetic data (NAD27 and NAD83) into the same grass location?
v.proj seems to be the way to do it, but I tried and it didn't work as I
expected it to.

Slightly more detailed version: 

The USGS has released two separate sets of files for the area I'm playing
with (Tijeras, NM).  The older set contains the DEM file, and a whole slew of
DLGs for roads, rivers, contours, etc.  The only new one that's available for
free in SDTS format is the roads.  The older roads file doesn't contain some
significant trails, and is fairly outdated --- Tijeras is a
semi-rural/suburban area that has changed a lot since the old DLG was
created.  Most of the other features (hypsography, hydrography, and of course
the DEM :-) haven't changed much in the same time.  Some roads have been
re-routed, others have been extended, and as I said hiking trails and jeep
roads appear in the newer one but not the older.

Thing is, the older data all used the NAD27 datum, and the new uses NAD83.
Both sets, however, cover a quadrangle labeled with the same lat/lon
coordinates (SW=35d00'00"N,106d30'00"W,  NW=35d7'30"N,106d30'00"W,
NE=35d7'30"N,106d22'30"W, SE=35d00'00"N,106d22'30"W).   The USGS 7.5' quad
map I have uses a UTM projection with the NAD83 datum, and that map has corner
ticks to represent where the same lat/lons would be placed with NAD27.
Eyeballing the scale says that these NAD27 ticks are about 50m to the left of
the NAD83 corners, with no visible difference in the vertical direction.

Ok, so what's the problem?  Grass has v.proj, and that should be able to
import my NAD83 DLGs from one location to another.  Well, I tried that.

When I create a grass location for the two different sets of data and use
v.proj to pull the roads map from the NAD83 location into the NAD27 location,
the result is that the outlines (lines of constant lattitude and longitude)
area mapped right onto one another, but the roads themselves are shifted by an
amount just about equal to the shift shown by the NAD27 corner ticks on my
quad map!  That is, if I go into the NAD27 location, display the roads vector
map, v.proj the NAD83 roads map in and redisplay, the two vector maps have the
same lat/lon outline, but none of the roads coincide.

I've tried to use the NADCON software from USGS, and I've convinced myself
that the lat/lon corners shouldn't map onto each other.  I'm now very
confused.  Has *anyone* seen this sort of thing before, and how have you fixed
it?  

[Longwinded details follow for those who want it.  You have been warned.]
Here's what I did:

I created two locations, tijeras27 and tijeras83.  I created each one by
selecting UTM projection and the appropriate ellipsoid and datum for the two
data sets.  Here's what the PROJ_INFO and DEFAULT_WIND files look like for the
two: 

tijeras27:
PROJ_INFO                   DEFAULT_WIND             
name: UTM                   proj:       1            
datum: nad27		    zone:       13           
dx: -22.000000		    north:      3887736.75   
dy: 157.000000		    south:      3873709.25   
dz: 176.000000		    east:       374713.062   
proj: utm		    west:       363113.9687  
ellps: clark66		    cols:       387          
a: 6378206.4000000004	    rows:       468          
es: 0.0067686580	    e-w resol:  30           
f: 294.9786982000	    n-s resol:  30           
zone: 13

I generated the tijeras27 location by the following method:
  - I started Grass by choosing the Spearfish location.
  - I used r.in.gdal to import the SDTS DEM file for Tijeras, NM and create a
    location.
  - I exited Grass and re-entered, selecting the new location.
  - the DEM imported by r.in.gdal is full of junk --- all the null values 
    between the edge of the rectangular UTM block and the (rotated) lat/lon
    quadrangle are set to -32766 or some such.  So I used r.out.ascii to 
    create an ascii file, used sed to change the -32766 to 0.  
  - Since r.in.gdal didn't create PROJ_INFO and PROJ_UNITS files when it
    created the tijeras27, I went ahead and exited grass5, deleted the
    tijeras27 directory, then created it over again through the normal "create
    new location" procedure, using the N,S,E and W corners in the header of
    the new ascii file.  I chose projection=utm, zone=13, ellipsoid=clark66
    datum=nad27 when I created the new location.
  - I re-imported the DEM using r.in.ascii
  - I imported all the NAD27 DLG files using v.in.sdts

Heres the relevant files from the other location, tijeras83:
PROJ_INFO                  DEFAULT_WIND           
name: UTM		   proj:       1          
datum: nad83		   zone:       13         
dx: 0.000000		   north:      3887934.73 
dy: 0.000000		   south:      3873906.74 
dz: 0.000000		   east:       374715.98  
proj: utm		   west:       363117.12  
ellps: grs80		   cols:       387        
a: 6378137.0000000000	   rows:       468        
es: 0.0066943800	   e-w resol:  30
f: 298.2572221010	   n-s resol:  30
zone: 13

I got these corner locations by using m.ll2u on the four corner lat/lons and
choosing the appropriate min/max values.  I gave the lat/lon pairs above, and
just re-verified that m.ll2u gives those values.

I then imported the new NAD83-referenced DLG for roads into this location.

Lastly, I went into the tijeras27 location:
v.proj input=roads location=tijeras83 output=roads83

Displaying them both with d.vect shows the roads themselves offset to the left
by something like 50-100 meters.  The outlines line up exactly.

To double check that they really were created with different datums and that
projection really was needed, I tried importing the newer DLG into the
tijeras27 location directly with v.in.sdts --- that was far worse, nothing
lined up.


All of this data that I'm using is freely available from the EROS home page at
http://edcwww.cr.usgs.gov/doc/edchome/ndcdb/ndcdb.html
I'm using the 1:24K DEM and 1:24K DLG files for Tijeras, New Mexico if you
care to duplicate this.  To completely reproduce what I have here it would be
enough to download just the two roads files 1036217.RD.sdts.tar.gz and
1452479.RD.sdts.tar.gz.  I gave the DEFAULT_WIND files, so you could easily
create the location without much trouble.


I said above that I convinced myself that the corners shouldn't line up by
using the NADCON software from USGS (ftp://kai.er.usgs.gov).  Here's what I
tried:

#Map the SW corner of the NAD83 quad into NAD27 UTMs:
GRASS:~ > nad2nad -i 83,rev -o 27,utm=13 -r conus
35d0'0"n 106d30'W
363166.56       3873867.07

#Map the NAD27 UTM for that back to NAD83 UTMS
GRASS:~ > nad2nad -i 27,utm=13 -o 83,utm=13 -r conus
363166.56       3873867.07
363117.13       3874070.93

#Verify that the new NAD27 coordinates really map back to 35d0'0"n 106d30'W
GRASS:~ > nad2nad -i 27,utm=13 -o 83  -r conus
363166.56       3873867.07
106d30'W        35dN

#Now try mapping the NAD83 UTMs to NAD27 LAT/LON

GRASS:~ > nad2nad -i 83,utm=13 -o 27,rev  -r conus
363117.13       3874070.93
34d59'59.82"N   106d29'57.922"W


Now lookie that.  Rather than mapping to (35d,106d30') they are a little off,
less than a second in lattitude, but about 2 seconds in longitude.  Now, that
could be conversion roundoff error.  So let's see what (35d,106d30') in NAD27
maps to in NAD83 utms: 

#convert 35dN,106d30'W from NAD27 to NAD83 UTM: 
GRASS:~> nad2nad -o 83,utm=13 -i 27,rev -r conus 
35d0'0"n 106d30'W 
363064.53	3874077.26

Now putting it all together:
NAD83 (35d0'0"n,106d30'W) = (363117.13Easting, 3874070.93Northing) NAD83
NAD27 (35d0'0"n,106d30'W) = (363064.53Easting, 3874077.26Northing) NAD83

and continuing that line of nad2nads (I won't include any more) we get
NAD27 (35d0'0"n,106d30'W) = (363113.96Easting, 3873873.40Northing) NAD27
NAD83 (35d0'0"n,106d30'W) = (363166.56Easting, 3873867.07Northing) NAD27

[Note: I just reread my DEFAULT_WIND files before sending this --- see how the
 SW corner's easting matches the "west" lines for the appropriate locations?
 Whew.  Just a little sanity check.  Lord knows I need one right now.]

So if I were to properly project the NAD83 quadrangle onto the same display as
the NAD27 quadrangle, the southwest corner of the NAD83 should be about
53meters to the right of the southwest corner of the NAD27 quad, and about 6
meters (neglible) down.  This is what the paper quad map shows (to within an
eyeball-norm, remember what I wrote in the earlier part of the message?).

So why doesn't v.proj produce that?  How is it that v.proj can say that the
two quadrangles line up, and then leave all the interior details to be
misaligned?  What am I missing?

Now, if you've made it this far, and you can help me, please, please, please
drop me a line (or say it here, so it'll get archived for others to benefit
from).  I'm actually playing with this stuff so much I'm not getting enough
sleep, so it's time to go to bed.  And then, of course, I have to start
figuring out what the deal is with v.reclass and v.reclass.pg that they won't
let me reclassify roads by type...but that's for another night just before
midnight.

-- 
Tom Russo   KM5VY   QRPL #1592   K2#398     http://www.swcp.com/~russo/
Tijeras, NM DM64ux  SOC #236  ICQ#97201722  http://www.qsl.net/~km5vy/
 "What if the hokey pokey really is what it's all about?"


