[GRASSLIST:1868] Importing USGS SDTS DLG maps from differing datums? I'm confused. [very long]

GRASS GRASS user
,
1

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 :slight_smile: 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?"

2

On Thu, 24 May 2001, Tom Russo wrote:

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.

Tom,

I think you found an actual problem. I'm nearby, so I took a local
interest in your problem.

I reproduced the problem that you had, then tried an alternative
procedure. I converted both sets of road data from UTM coordinates to the
corresponding lat-long coordinates, then tried to convert the NAD83 roads
from lat-long NAD83 to lat-long NAD27 so the two data sets could be
compared.

This procedure parallels the process that you have to go through to
convert the data using the National Coast and Geodetic Survey's (NCGS)
software. The conversion from lat-long NAD83 to lat-long NAD27 failed.
v.proj produced:

Creating dig file ...
Error in pj_do_proj

and the output file was empty.

I believe you were correct that the boundaries of the two maps should
*not* have overlain after the conversion. So I did a complete conversion
manually to check GRASS' results.

I dumped the NAD83 roads to an ascii vector file and wrote a small program
that converted the data to a format that could be read by the NCGS
programs. I ran that data through th NCGS' "utms" program to convert from
UTMs to lat-long on the GRS80 spheroid, then used the NCGS NDCON210
program to convert the latitudes and longitudes from NAD83 to NAD27. I
again used "utms" to convert the NAD27 lat-long records back to UTM
coordinates on the Clark '66 spheroid. Finally, I wrote a small program
that merged the converted coordinates back into the original ascii vector
file.

I've been using this procedure to convert NAD83 UTM values to NAD27 UTM
values for several years, and I'm pretty sure that I got it right.

I imported the new ascii vector file back into GRASS and compared it to
the USGS' NAD27 roads. The two files were pretty much the same except for
the obvious differences caused by new road development.

So, I think you have found a fairly serious bug in GRASS. Do you want to
do the official bug report?

As a product of all that I have this ascii vector file that has the new
roads system in NAD27 UTMs. If you want it just send me an email.

Roger Miller
Lee Wilson and Associates

3

On Thu, May 24, 2001 at 03:38:52PM -0600, Roger S. Miller wrote:

On Thu, 24 May 2001, Tom Russo wrote:

> 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.

Tom,

I think you found an actual problem. I'm nearby, so I took a local
interest in your problem.

Shortcoming actually. Don't think GRASS actually can handle datum
transformations yet (though there's no reason it shouldn't). I believe
Andreas added the ability to specify a datum, and some functions that
will do datum transformations, but these haven't been tied into existing
routines for inverse and forward projections. I think do_proj.c
(libes/proj/) will need to be hacked to do datum transformations when the
input/output datums are different. It simply doesn't look, and libproj
does not do datum shifts.

--
Eric G. Miller <egm2@jps.net>

4

"Eric G. Miller" wrote:

Shortcoming actually. Don't think GRASS actually can handle datum
transformations yet (though there's no reason it shouldn't). I believe
Andreas added the ability to specify a datum, and some functions that
will do datum transformations, but these haven't been tied into existing
routines for inverse and forward projections. I think do_proj.c
(libes/proj/) will need to be hacked to do datum transformations when the
input/output datums are different. It simply doesn't look, and libproj
does not do datum shifts.

Thanks for the insight, Eric. I see you're right. When pj_do_proj is
asked to convert lat-long to lat-long it simply returns an error. It
doesn't check for the datum at all.

OK, I can accept that as a shortcoming in the code, but it's a pretty
big shortcoming when you consider what's going on in North America right
now.

Although it isn't necessarily a bug in the code, it's definitely a bug
in the documentation. The user is prompted for the datum when he
creates a location and g.region shows the datum as part of the location
definition. The documentation for v.proj says that the v.proj converts
between mapsets in different locations so I think most GRASS users will
logically assume that v.proj converts between different datums.

At a minimum the documentation needs to be changed. Probably the code
should be changed to detect when a user is trying to convert between
datums and return an appropriate statement.

I have to wonder how many people have already used GRASS for that
conversion and now have erroneous data as a result. I'm pretty sure
that I've done it.

Roger Miller

5

On Thu, May 24, 2001 at 09:50:22PM -0600, Roger S. Miller wrote:

Thanks for the insight, Eric. I see you're right. When pj_do_proj is
asked to convert lat-long to lat-long it simply returns an error. It
doesn't check for the datum at all.

OK, I can accept that as a shortcoming in the code, but it's a pretty
big shortcoming when you consider what's going on in North America right
now.

Although it isn't necessarily a bug in the code, it's definitely a bug
in the documentation. The user is prompted for the datum when he
creates a location and g.region shows the datum as part of the location
definition. The documentation for v.proj says that the v.proj converts
between mapsets in different locations so I think most GRASS users will
logically assume that v.proj converts between different datums.

At a minimum the documentation needs to be changed. Probably the code
should be changed to detect when a user is trying to convert between
datums and return an appropriate statement.

I have to wonder how many people have already used GRASS for that
conversion and now have erroneous data as a result. I'm pretty sure
that I've done it.

I don't know how many, but Andreas only completed these changes late
last fall / early winter. Definitely a big sign in the documentation
should say that datum transformations aren't handled. Better to hack
the projection routines, though I think there we'll have to make sure
"coorconv" gets built before "proj". I don't know how to determine
which method should be used for datum transformation (usually
3-parameter, I guess). We still don't carry NADCON files, so nad27 to
nad83 will be less than ideal anyway. The "proj" distribution comes with
a nad2nad utility, but it doesn't seem to be integrated. I've packaged
up a library from NIMA's GeoTrans calculator which would handle datum
shifts rather transparently. However, it doesn't use NADCON files
either and I think the number of supported projections is more limiting
(25, if I remember).

A while ago, a problem with state plane was brought up, and Morten and I
discussed some major reorgs of the g.setproj as well as modifying some
of the support routines (for instance, a datum implies an ellipsoid so
the user shouldn't be bothered with both). We agreed to hold off for
5.1 because there've been enough hold ups to getting 5.0 out.

I would like to see GRASS adopt the EPSG system for identifying
projections, datums, ellipsoid (and possibly entire coordinate systems)
as this is used by GeoTIFF and a number of commercial vendors are also
migrating to supporting this system. Since each of those things has a
numeric identifier it should be easy to do lookups quickly when needed.

--
Eric G. Miller <egm2@jps.net>

6

On Thu, May 24, 2001 at 09:58:43PM -0700, Eric G. Miller wrote:

On Thu, May 24, 2001 at 09:50:22PM -0600, Roger S. Miller wrote:
> Thanks for the insight, Eric. I see you're right. When pj_do_proj is
> asked to convert lat-long to lat-long it simply returns an error. It
> doesn't check for the datum at all.
>
> OK, I can accept that as a shortcoming in the code, but it's a pretty
> big shortcoming when you consider what's going on in North America right
> now.
>
> Although it isn't necessarily a bug in the code, it's definitely a bug
> in the documentation. The user is prompted for the datum when he
> creates a location and g.region shows the datum as part of the location
> definition. The documentation for v.proj says that the v.proj converts
> between mapsets in different locations so I think most GRASS users will
> logically assume that v.proj converts between different datums.
>
> At a minimum the documentation needs to be changed. Probably the code
> should be changed to detect when a user is trying to convert between
> datums and return an appropriate statement.
>
> I have to wonder how many people have already used GRASS for that
> conversion and now have erroneous data as a result. I'm pretty sure
> that I've done it.

I don't know how many, but Andreas only completed these changes late
last fall / early winter. Definitely a big sign in the documentation
should say that datum transformations aren't handled. Better to hack
the projection routines, though I think there we'll have to make sure
"coorconv" gets built before "proj". I don't know how to determine
which method should be used for datum transformation (usually
3-parameter, I guess). We still don't carry NADCON files, so nad27 to
nad83 will be less than ideal anyway. The "proj" distribution comes with
a nad2nad utility, but it doesn't seem to be integrated. I've packaged
up a library from NIMA's GeoTrans calculator which would handle datum
shifts rather transparently. However, it doesn't use NADCON files
either and I think the number of supported projections is more limiting
(25, if I remember).

A while ago, a problem with state plane was brought up, and Morten and I
discussed some major reorgs of the g.setproj as well as modifying some
of the support routines (for instance, a datum implies an ellipsoid so
the user shouldn't be bothered with both). We agreed to hold off for
5.1 because there've been enough hold ups to getting 5.0 out.

I would like to see GRASS adopt the EPSG system for identifying
projections, datums, ellipsoid (and possibly entire coordinate systems)
as this is used by GeoTIFF and a number of commercial vendors are also
migrating to supporting this system. Since each of those things has a
numeric identifier it should be easy to do lookups quickly when needed.

Hi,

a short remark on the GRASS 5.0.0 and no (not yet) datum support:

At time GRASS 5.0.0 is using PROJ 4.3 software from USGS without datum
transformation support. I agree that the query for datum is confusing,
at least a hint needs to be added that this feature will be used from 5.1
onwards.

From GRASS 5.1 onwards GRASS will use the latest PROJ 4.x software available

from http://www.remotesensing.org/proj/ maintained by Frank Warmerdam
including datum transform. As far as I know the "coorconv" library won't ne
necessary then as PROJ4.4.x already supports datum transforms (from the web
page: "Support has also been added for 3 and 7 parameter datum shifts").

Again, this is forthcoming for 5.1.

Regards

Markus Neteler

7

On Thu, 24 May 2001, Roger S. Miller wrote:

I believe you were correct that the boundaries of the two maps should
*not* have overlain after the conversion. So I did a complete conversion
manually to check GRASS' results.

Roger,

  You have what appears to be a complete description of the conversion
algorithm in your message. Would you be willing to make your untility
programs available to the rest of us who have to cope with mixed-data
(NAD27, NAD83), too?

  Many areas use both older and newer data. In my former life, working with
MapInfo, I converted everything to NAD83/feet (because that's what many
regulatory agencies want) and the applicable State Plane Coordinate zone.

Thanks,

Rich

Dr. Richard B. Shepard, President

                       Applied Ecosystem Services, Inc. (TM)
            2404 SW 22nd Street | Troutdale, OR 97060-1247 | U.S.A.
+ 1 503-667-4517 (voice) | + 1 503-667-8863 (fax) | rshepard@appl-ecosys.com