[GRASSLIST:1995] Re: v.proj bug?

[Paul Kelly]

Hello

David Orme wrote:
> 
> Hi,
> 
> I e-mailed a while back about a problem with s.proj - I've encountered
> something similar in v.proj.
> 
> I have a raster dataset in the old New Zealand South Island grid
> (Transverse Mercator: central meridian 171.5E, latitude of origin 44S,
> scale 1, false easting 500000, false northing 500000, units yards,
> ellipse international). This projection is included in ArcMAP and when
> I compare a graticuled map in ArcMAP to (coarse scaled) printed grid of
> the raster  for the dataset, there is a very close match.
> 
> I've tried to set up a GRASS location to match this (I've appended the
> PROJ_INFO) and downloaded the DCW vector coastline to use as a base
> map. I imported this into a latlong location and v.proj'd it across and
> found that the coordinates did not match ArcMAP. I've picked two South
> Island locations as test sets  - Cape Providence (bottom left) and the
> tip of Farewell Spit (top right). These have the following lat longs
> (in s.in.ascii format and cs2cs):
> 
> 166:27:28.0E 46:00:55.0S Cape_Providence
> 173:01:05.1E 40:33:28.8S Farewell_Spit
> 
> 166d27'28.0E 46d00'55.0S
> 173d01'05.1E 40d33'28.8S
> 
> In ArcMAP (and through a separate NZ grid application for Classic Mac
> and approximately from my printed grid) these have the NZ SI
> coordinates:
> 
> 72976E 241337N Cape_Providence
> 640603E 916716N Farewell_Spit
> 
> I've run the vector file and these two coordinates through cs2cs: the
> vector file imported through cs2cs and the e00 imported through v.proj
> overlie exactly and the coordinates of the points out of cs2cs are:
> 
> 119772.29E 288333.97N Cape_Providence
> 687425.44E 963732.29N Farewell_Spit
> 
> The cs2cs arguments used were:
> 
> cs2cs +proj=latlong +ellps=clrk66 +to +proj=tmerc +lat_0=44S
> +lon_0=171.5E +x_0=500000 +y_0=500000 +ellps=intl +units=yd

The false easting and northing should be in metres even though the
projected unit is yards; at 0.9144 yards to the metre this gives

cs2cs +proj=latlong +ellps=clrk66 +to +proj=tmerc +lat_0=44S \
+lon_0=171.5E +x_0=457200 +y_0=457200 +ellps=intl +units=yd

which comes out at
72965.64        241527.32
640618.79       916925.64
i.e. a lot closer to the values you were looking at.

The other issue (which should remove the remaining error when addressed)
is that you are projecting onto a different ellipsoid so you will need
to do a datum transformation. GRASS (version >=5.3) now supports the
NZGD49 datum (on which your map appears to be based) quite well; try the
following for your PROJ_INFO:

name: Transverse Mercator
proj: tmerc
datum: nzgd49
nadgrids: nzgd2kgrid0005.gsb
ellps: international
a: 6378388.0000000000
es: 0.0067226700
f: 297.0000000000
lat_0: -44.0000000000
lon_0: 171.5000000000
k_0: 1.0000000000
x_0: 457200.0
y_0: 457200.0

However it is not clear what datum your 'DCW' data is in but you will
also need to provide transformation parameters for that location. These
parameters should be suitable for transforming from the clark66
ellipsoid to the GRS80/WGS84 over your area of interest. When
re-projecting then GRASS/PROJ will go to the WGS84 as an intermediate
stage between clark66--->international.

Can you find what datum and parameters are being used for the clark66
data, from ArcMAP or your NZ grid Mac application?

I hope this isn't too late and you get something working.

Paul K