[GRASS5] GRASS, PROJ and NAD

neteler · October 11, 2003, 8:55am

I was playing a bit with PROJ and have generated
a Spearfish location with UTM/NAD83. Reprojection from
the sample NAD27 works well. Then I have created a
Spearfish-Lat-Long location and tried to reproject
data from the sample Spearfish location into.

Here I get the error:
r.proj in=elevation.dem location=spearfish mapset=PERMANENT
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
[...]

Being no NAD expert, where is the trick?

Thanks in advance,

Markus

Paul_Kelly · October 11, 2003, 1:16pm

On Sat, 11 Oct 2003, Markus Neteler wrote:

I was playing a bit with PROJ and have generated
a Spearfish location with UTM/NAD83. Reprojection from
the sample NAD27 works well. Then I have created a
Spearfish-Lat-Long location and tried to reproject
data from the sample Spearfish location into.

Here I get the error:
r.proj in=elevation.dem location=spearfish mapset=PERMANENT
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
[...]

That message usually means the point you are trying to re-project falls
outside the region covered by the NAD correction file, but there can be
several other reasons. Did r.proj not print the parameters for each
location first? What are the PROJ_INFO files?

You could experiment with s.proj just projecting one site at a time and
see if the co-ordinates made sense.

neteler · October 13, 2003, 9:56am

On Sat, Oct 11, 2003 at 02:16:28PM +0100, Paul Kelly wrote:

On Sat, 11 Oct 2003, Markus Neteler wrote:

> I was playing a bit with PROJ and have generated
> a Spearfish location with UTM/NAD83. Reprojection from
> the sample NAD27 works well. Then I have created a
> Spearfish-Lat-Long location and tried to reproject
> data from the sample Spearfish location into.
>
> Here I get the error:
> r.proj in=elevation.dem location=spearfish mapset=PERMANENT
> pj_transform() failed
> cause: failed to load NAD27-83 correction file
> pj_transform() failed
> cause: failed to load NAD27-83 correction file
> pj_transform() failed
> cause: failed to load NAD27-83 correction file
> pj_transform() failed
> cause: failed to load NAD27-83 correction file
> pj_transform() failed
> [...]

That message usually means the point you are trying to re-project falls
outside the region covered by the NAD correction file, but there can be
several other reasons. Did r.proj not print the parameters for each
location first? What are the PROJ_INFO files?

Ah, sorry:

In spearfishLL:
r.proj in=elevation.dem location=spearfish mapset=PERMANENT

Input Projection Parameters: +proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +nadgrids=conus
Input Unit Factor: 1

Output Projection Parameters: +proj=latlong +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000
Output Unit Factor: 1
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
pj_transform() failed

[ deleted 450000 lines ]

cause: failed to load NAD27-83 correction file
pj_transform() failed
cause: failed to load NAD27-83 correction file
Input:
Cols: 447 (633)
Rows: 104 (466)
North: 4920620.000000 (4928000.000000)
South: 4917500.000000 (4914020.000000)
West: 595590.000000 (590010.000000)
East: 609000.000000 (609000.000000)
ew-res: 30.000000
ns-res: 30.000000

Output:
Cols: 61 (129582)
Rows: 10 (10)
North: 44.433333 (44.433333)
South: 44.405556 (44.405556)
West: -103.800000 (-103.800000)
East: -103.630556 (256.150000)
^^^^^^^^^^^^^^^^^
ooops!! I just realize it now.
ew-res: 0.002778
ns-res: 0.002778
Allocating memory and reading input map...100%
Projecting... 100%

There is something strange. Here is the spearfishLL PROJ_INFO file:

name: Latitude-Longitude
datum: wgs84
towgs84: 0.000,0.000,0.000
proj: ll
ellps: wgs84

Here is the spearfishLL DEFAULT_WIND file:
proj: 3
zone: 0
north: 44:26N
south: 44:24:20N
east: 103:51W
west: 103:48W
cols: 129582 <--??
rows: 10
e-w resol: 0:00:10
n-s resol: 0:00:10

The cols value is surprising (to me). However, this is what the
dialogue provides when defining spearfishLL:

  projection: 3 (Latitude-Longitude)
  zone: 0
  north: 44:26N
  south: 44:24:20N
  east: 103:51W
  west: 103:48W

e-w res: 0:00:10
n-s res: 0:00:10

  total rows: 10
  total cols: 129582
  total cells: 1,295,820

Do you accept this region? (y/n) [y] >

Maybe I am missing something, but above looks unusual.

Markus

H_B · October 13, 2003, 10:28am

Here is the spearfishLL DEFAULT_WIND file:
proj: 3
zone: 0
north: 44:26N
south: 44:24:20N
east: 103:51W
west: 103:48W
cols: 129582 <--??
rows: 10
e-w resol: 0:00:10
n-s resol: 0:00:10

The cols value is surprising (to me). However, this is what the
dialogue provides when defining spearfishLL:

The 'east' is actually west of 'west', so that the longitudinal extent wraps around the planet, for all but the 5km between the two!

129582 cols * 1/6 nmile = 21597.

$ gcb n 44.25N103.51W 44.25N103.44W
Short Path Bearing is 090 Degrees for 5 Nautical Miles
Long Path Bearing is 270 Degrees for 21595 Nautical Miles

.. swap west and east ..

Hamish

H_B · October 13, 2003, 10:52am

The 'east' is actually west of 'west', so that the longitudinal extent
wraps around the planet, for all but the 5km between the two!

129582 cols * 1/6 nmile = 21597.

$ gcb n 44.25N103.51W 44.25N103.44W
Short Path Bearing is 090 Degrees for 5 Nautical Miles
Long Path Bearing is 270 Degrees for 21595 Nautical Miles

and that's a bug in gcb.

R(lat) = Rearth*cos(90-lat)

R(44.417) = 6367.4km * cos(90 - 44.417) = 4456.3km

Circumference along the 44.417 parallel is therefore 4456.3km * 2PI
= 28,000km or 15,119 nautical miles.

As to what effect that will have in GRASS on global datasets, such as
the accuracy of data viewed with d.*.viewproj ....? Not much we can do
about it I guess without going to a 3d grid.

Hamish

H_B · October 13, 2003, 10:56am

> 129582 cols * 1/6 nmile = 21597.

...

As to what effect that will have in GRASS on global datasets, such as
the accuracy of data viewed with d.*.viewproj ....? Not much we can do
about it I guess without going to a 3d grid.

never mind. not thinking. 1/6 nmile should be 1/6' and is therefore scaled..

great circle calculator:
http://packages.debian.org/gcb

neteler · October 13, 2003, 11:47am

On Mon, Oct 13, 2003 at 11:28:10PM +1300, Hamish wrote:

>
> Here is the spearfishLL DEFAULT_WIND file:
> proj: 3
> zone: 0
> north: 44:26N
> south: 44:24:20N
> east: 103:51W
> west: 103:48W
> cols: 129582 <--??
> rows: 10
> e-w resol: 0:00:10
> n-s resol: 0:00:10
>
> The cols value is surprising (to me). However, this is what the
> dialogue provides when defining spearfishLL:

The 'east' is actually west of 'west', so that the longitudinal extent
wraps around the planet, for all but the 5km between the two!

Ah, good point!
Now it works out of the box

Wasn't there once a test to catch that errors?

[...]

Thanks for the hint, now r.proj is happy as well.

Markus

neteler · October 13, 2003, 11:57am

On Mon, Oct 13, 2003 at 01:47:39PM +0200, Markus Neteler wrote:

On Mon, Oct 13, 2003 at 11:28:10PM +1300, Hamish wrote:
Wasn't there once a test to catch that errors (W vs E)?

[...]
Another comment:
pj_transform() failed
cause: failed to load NAD27-83 correction file
This message (from ./proj/pj_strerrno.c?) is not really
helpful, should it be changed to "out-of-region" error?

Markus