From: gie@charon.er.usgs.gov (Gerald I. Evenden)
Newsgroups: comp.infosystems.gis
Subject: How do you gctp and proj to "do thier thing"???
Date: 15 Sep 92 18:51:31 GMT
Organization: U.S. Geological Survey

> >From:         Steve Murphy <murf@MURPHY.SPS.MOT.COM>
> Subject:      How do you gctp and proj to "do thier thing"???
> 
> Situation: For fun, I'm trying to play with the usgs 1:2M data off the
> CD  Rom. (You can pick it  off spectrum.xerox.com, if  you have the $$
> for  the disk and  net bandwidth). Anyway,  Before  I can really start
> having   fun, I need to  convert   the &^%%$^ albers-equal-Area coords
> internal to the file to lat-long data.
> 
> I used f2c on gctp-II that I got from usgs. No  dice. Doesn't run, and
> even when  I tweak, it  doesn't seem  to work  right.  Bad numbers all
> over. I call out the  initialization and forward and reverse  routines
> in pj03xxx.c (not looking at the code right now). I base my statements
> on the stuff pumped out to the error file.

The version of f2c that I employ compiles both the library and sample
programs with no problems.  Sample runs will not perform because of special
flags they used on NMD's host system which maps blanks internal to a
D field to 0's.  In standard FORTRAN they should be skipped.  Unless you
have experience with FORTRAN and the FORTRAN programmer's mentality,
this stuff is confusing.

National Geodetic Survey's version of GCTP is somewhat cleaner but no
driver programs are supplied and it does not have some of the newer
projections recently added to the USGS version.

> I resorted to proj,  from usgs. My  mapset has this header (split into
> 80chars/line:)
> 
> USGS-NMD  DLG DATA - CHARACTER FORMAT - 09-29-82 VERSION
> NORTHERN GREAT LAKES STATES              1967, 1972  2000000.
> 
>      3     3  9999     2 0.50800000000D+02     4     0     4     1
>    0.637820640000000D+07   0.676865799729109D-02   0.290300000000000D+08
>    0.450300000000000D+08  -0.960000000000000D+08   0.230000000000000D+08
>    0.0                     0.0                     0.0
>    0.0                     0.0                     0.0
>    0.0                     0.0                     0.0
>  0.10000000000D+01 0.0               0.0               0.0
> SW       42.000000  -97.000000         -82254.11  2111326.60
> NW       48.000000  -97.000000         -75202.08  2778443.33
> NE       48.000000  -82.000000        1049619.17  2855498.19
> SE       42.000000  -82.000000        1147578.85  2195575.19

I believe the x-y values are based upon the digitized corners of the
map and are only good to .001" (1:2,000,000 scale, see below).
The true precision of these numbers is not as published.

> POLITICAL BOUNDARIES   0   560   560 010   247   247 010   800   800   1
> N    1   -82254.11  2111326.60           0           0     0
> N    2   -75202.08  2778443.33           0           0     0
> N    3  1049619.17  2855498.19           0           0     0
> N    4  1147578.85  2195575.19           0           0     0
> N    5   -37075.92  2278381.19           2           0     0
> 
> Proj ALMOST works. The numbers are close,  but definitely OFF. Look at
> the lat/long vs local coords with the SW/NW/NE/SE entries.
> 
> Projs says:
> 
> ? proj +proj=aea +lon_0=96w +lat_0=23 +ellps=clrk66 +lat_1=29.5 +lat_2=45.5
> -97 42
> -82267.02	2111279.30
> -82 48
> 1049656.73	2855547.18

Looks good to me and I get the same numbers here.

> It appears that clrk66, lat_1=29.5, and lat_2=45.5 are all defaults in
> proj. Most of the numbers that gctp map into args for proj:
> 
>    0.637820640000000D+07   0.676865799729109D-02   0.290300000000000D+08
>    0.450300000000000D+08  -0.960000000000000D+08   0.230000000000000D+08
>    0.0                     0.0                     0.0
> 
> the first is the ellps stuff in proj:
>      GRS67: a=6378160.0      rf=247.247167    GRS 67(IUGG 1967)
>     bessel: a=6377397.155    rf=299.152813    Bessel 1841
>     clrk66: a=6378206.4      b=6356583.8      Clarke 1866
>     clrk80: a=6378249.145    rf=293.4663      Clarke 1880 mod.
>    everest: a=6377276.3452   b=6356075.4133   Everest 1830

Note: the proj +ellps list is parameters are based upon what appears
to be the defining numbers for the ellips.  In the case of Clarke 1866
the major and minor semi-radia are usually specified.  Many of the other
ellipsoid use a and 1/r.

> the a=6378206.4 lines right up with the data[0] in the list above.
> the second  number ? "Eccentricity^2 of  ellipsoid" this and #1 should
>                       line up with Clarke 1866.
> the 29d30m0s lines up with the lat_1 parameter of proj
> the 45d30m0s lines up with the lat_2 parameter of proj
> the -96d lines up with the lon_0 parameter
> the 23d lines up with lat_0 parameter.
> 
> The numbers come out close:
> -97 = -82267.02 in proj and usgs dataset said -82254.11
> -82 = 1049656.73 in proj and usgs dataset said 1049619.17

Your selection of proj parameters is correct.  The differences you
get between the ROM data and proj all appear to be less than 50m 
(see below).

> Can I do any better on accuracy with this? I was kind of hoping to get
> within +-.01 or so, not  +-50!!  It says in  the  data guide that  the
> resolution is 50.8 meters, and I'm not sure  how big this would  be in
> local coordinates,  but I  doubt much  left   of the  decimal point...
> Besides,  this resolution applies   to  actual  positional data in the
> field....   Why   would  I    be  off at     all  on   the  coordinate
> transformations?

The data for this CD-ROM was obtained by digitizing 1:2,000,000 scale
maps (with, I believe, a scanner).  Resolution accuracy of this process
was 0.001" or, as the file states, 50.8m for 1:1 scale.  I suspect that
the real accuracy is somewhat lower, especially when factoring in the
original accuracy of the map base scanned.

I would not recommend using the data for any larger scale efforts.

> Any proj experts or  people who've played with  gctp like to  speak up
> and tell me  how big a fool I  am? I  could stand to  learn something,
> that's for sure!

For better or worse, I am the author of program proj---next best thing
to finding an expert.

You are only misled by the common practice of tacking more significant
figures onto a number than are truly justified and believing everything
you read.  Naive, maybe.
Or as friend used to say: "it's not that I didn't know, I never suspected."

You did fine, given the information you had.

> murf

My door's always open.

G.I.Evenden   gie@charon.er.usgs.gov
Comments: Gated by NETNEWS@AUVM.AMERICAN.EDU