Newsgroups: comp.infosystems.gis
Date:         Mon, 14 Sep 1992 11:21:35 CDT
              <GIS-L@UBVM.BITNET>
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.

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

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

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

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?

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!

murf
Comments: Gated by NETNEWS@AUVM.AMERICAN.EDU