[GRASSLIST:1618] Re: More projection problems...

[Morten Hulden]

Hi Matt (and others following this thread),

I have compared the projection results given by GRASS' internal PROJ
with results given by external PROJ (4.3.3) using your numbers. They
match. So GRASS in projecting exactly as external PROJ both for aea->ll
and stp->ll. 

But I have some doubts here about the accuracy of the aea map parameters:

Where is the central parallel in your aea map? 

Your earlier message mentions 24N (same as the 1st standard parallel), but
in your PROJ_INFO file you have it at 27.45 (Did you mean to put it at
27.75, which is exactly midway between the two standard parallels, and
wrote 27.45 instead of 27:45?) Anyway, the latitude of the central
parallel affects the aea->ll projektion a lot. With the central parallel
at 24N you get the latitude of your test point at roughly the expected
29N, but with 27.45 you get it at 32N. Way outside the borders of your STP
location. r.proj is correct when it claims the input area is outside the
output region!

How accurate are the parameters given for the aea map?

Try using a false easting of 355000 instead of 400000 and your test points
longitude will get very close to the expected 81:53W. 

It also puzzles me that you have an Albers Equal Area map with a one-meter
resolution. Albers is mostly used for displaying larger areas, like whole
continents, with at most 500-1000 m resolution. In 1m resolution the
accuracy of the parameters start affecting the result a lot - the accuracy
should be roughly at the same level as the resolution. And your input
location is only about 7000x7000 feet in size with a 1 foot resolution.
Aren't you expecting too much from the Albers map?

I note that you use the upper left corner of both locations as the test
point and that you expect them to match. What is your reason for that
assumption? I would try to identify some features with known
co-ordinates on the map and use them as test points. Not just one point -
but several. As Eric mentioned, i.rectify is another tool in Grass that
may be more suited for the task you are trying to perform, especially when
the parameters of the input map are unknown or inaccurate.

Finally, since Albers is a conical and STP is a cylindrical projection,
you will not get a rectangular area when you project maps between the
locations. The projected maps will be rotated and their shapes
distorted, and the projected input may not cover the output location
completely if you tried to match the locations by a corner point. Better
try to define the locations so the _centers_ of the maps are
roughly in the same geographic area (and make the input a little larger
than the output location).

Tip: in order to speed up test projections you can define regions with
lower than the default resolution within your output location. r.proj is
region sensitive for resolutions too!. Then, when you are ready for real
work, select your default region and reproject the map.


best regards
Morten Hulden


On Thu, 15 Mar 2001, mberglund wrote:

> > I get exactly the same numbers 380339.35 1764858.00 using the same
> > parameters you have (+ ellipsoid=wgs84)
> I have been using wgs84 also.
> 
> > You can use m.proj to break down the projection into two steps:
> Ok, I'm with you. Lets do a couple of steps together.
> 
> First the X and Y below are the ACCURATELY CONVERTED STP coordinates for
> the original point in aea. I did the conversion to ll from his numbers. So
> these are the expected coordinates for the other rounded point.
> 
>  X (feet)        Y (feet)        Longtitude      Latitude
> -------------   --------------  -----------     ---------
> 
>  536507.00      1764928.00    81:53:08.131239W 29:11:20.067269N
> 
> By the way, my original estimate was wrong, I'm off by more like 156000
> feet. Yikes!!
> 
> > My m.proj gives 
> > 82d22'30.036W 29d11'17.659N 
> > give (i assumed the ellipsoid wgs84). Does that make sense?
> Ok. I got the ellipsoid(s) he is using (he wasn't certain of which but it
> is one of the two) grs80 or wgs84. Stephen at FSU works on a number of
> projects, so he couldn't remember which one off the top of his head. 
> I've tried both, so I doubt that this is an issue.
> 
> > Note that r.proj nor m.proj will not make datum.shifts (yet), but datum
> > differences should not give errors of the magnitude you mention.
> Fortunately I now we know we don't have to do datum shifts.
> 
> > Also note that the proj parameters for STP projections are statically
> > defined in <grass install directory>/etc/state27 and FIPS.code files. I
> > have no idea about the accuracy of those formulas. Maybe someone with more
> > knowledge about STP would like to check them.
> I am not very knowledgeable, but as I understand it, these numbers are
> gotten from PROJ.4. Those numbers have been good enough for me to go from
> lat/long to stp from gps data, and present it on a map within a couple
> foot accuracy. 
>  
> > Can you dump us your PROJ_INFO for both the input and output location, the
> > DEFAULT_WIND for the output location and the cellhd file of the input map
> > (a file with the same name as the input map, located in the cellhd
> > directory of the input location.) Also please say what your expected STP
> > co-ordinates for your test poit are exactly. This would help to do some
> > more debugging.
> 
> The test point expected results(also in the chart above) are:
> X: 	536507.00
> Y:	1764928.00
> NOTE: The original points were rounded for expediance. The resultant
> points that were generated from the rounded points were still thousands of
> feet more accurate than my conversion. 
> 
> STP Location:-------------------------
> DEAFULT_WIND-
> 	proj:       99
> 	zone:       0
> 	north:      1907628
> 	south:      1657863
> 	east:       784820
> 	west:       473816
> 	cols:       311004
> 	rows:       249765
> 	e-w resol:  1
> 	n-s resol:  1
> 
> PROJ_INFO-
> 	name: State Plane
> 	proj: tmerc
> 	a: 0.63782064e+07
> 	es: 0.6768657997291094e-02
> 	x_0: 0.1524003048006096e+06
> 	y_0: 0
> 	k: 0.9999411764705882e+00
> 	lon_0: 82dw
> 	lat_0: 24d20'n
> 
> Aea Location:---------------------------
> DEFAULT_WIND-
> 	proj:       99
> 	zone:       0
> 	north:      576795
> 	south:      569695
> 	east:       564255
> 	west:       557754
> 	cols:       6501
> 	rows:       7100
> 	e-w resol:  1
> 	n-s resol:  1
> 
> PROJ_INFO-
> 	name: Albers Equal Area
> 	proj: aea
> 	ellps: wgs84
> 	a: 6378137.0000000000
> 	es: 0.0066943800
> 	f: 298.2572235630
> 	lat_0: 27.4500000000
> 	lat_1: 24.0000000000
> 	lat_2: 31.5000000000
> 	lon_0: -84.0000000000
> 	x_0: 400000.0000000000
> 	y_0: 0.0000000000
> 
> Hope this helps and Thanks again,
> Matt
>