[GRASS-user] scale error for each pixel

[Ken Mankoff]

Dear All,

Thank you for the suggestions. In case this is useful for anyone else, I'll give my (current) solution below.

On 2017-11-21 at 08:30, Markus Metz <markus.metz.giswork at gmail.com> wrote:
> On Tue, Nov 21, 2017 at 9:26 AM, Moritz Lennert <
> mlennert at club.worldonline.be> wrote:
>
> You could run area() both in your projected location and in a lat-long
> location (preferably using a similar datum), although you would have
> to be careful in the definition of the resolution in each location.

A good idea but I'd be wary, as you point out, of getting the resolution correct. This may be complicated by my region covers large areas near the poles.

> You can avoid problems with raster resolution by creating a vector
> grid with areas. Calculate the area sizes for the vector grid, then
> reproject the vector grid and calculate the area sizes again. Maybe
> the differences in area sizes for the same vector area give you what
> you are looking for.

Yes. But I think you suggested checking with the proj mailing list and they gave good advice that seems easier than the above suggestion *for small areas*. I think if I want a raster grid with errors everywhere, your suggestion above might be the way to go. Or perhaps just pass a few 100s of 1000s of points to the proj or geod command below and then read ASCII file of errors directly in to a raster.


1) proj command with "-VS" flag gives lots of info:

# test mapset
grass72 -c EPSG:3413 -e ./Gproj
grass72 ./Gproj/PERMANENT

# Greenland
g.region  n=-632500 s=-3384500 w=-653000 e=879800 res=100 -a
g.region -p

PROJSTR=$(g.proj -j)
echo $PROJSTR

GRASS 7.2.0 (Gproj):> echo "-45 60" | proj -VS ${PROJSTR}
#Stereographic
#	Azi, Sph&Ell
#	lat_ts=
# +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +no_defs
# +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +to_meter=1
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257224
#  squared eccentricity: 0.006694379990
Longitude: 45dW [ -45 ]
Latitude:  60dN [ 60 ]
Easting (x):   0.00
Northing (y):  -3323160.27
Meridian scale (h) : 1.03942808  ( 3.943 % error )
Parallel scale (k) : 1.03942808  ( 3.943 % error )
Areal scale (s):     1.08041073  ( 8.041 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ -0.00000000 ]
Max-min (Tissot axis a-b) scale error: 1.03943 1.03943



I assume because the above is using the GRASS proj string that the error reported here really are map projection errors, and not due to some datum or geoid shift.



2) Geod command

# test mapset
grass72 -c EPSG:3413 -e ./Gproj
grass72 ./Gproj/PERMANENT

# Greenland
g.region  n=-632500 s=-3384500 w=-653000 e=879800 res=100 -a
g.region -p

ew0=$(echo "879000|-2632000"  | m.proj -o input=-)
echo $ew0
e0=$(echo $ew0 | cut -d'"' -f1)
n0=$(echo $ew0  | cut -d"|" -f2 | cut -d'"' -f1)

# shift by 1000 m
ew1=$(echo "878000|-2632000"  | m.proj -o input=-)
echo $ew1
e1=$(echo $ew1 | cut -d'"' -f1)
n1=$(echo $ew1  | cut -d"|" -f2 | cut -d'"' -f1)

geod +ellps=WGS84 << EOF -I +units=m
  ${n0}N ${e0}W ${n1}N ${e1}W
EOF

>> -71d32'0.019"	108d26'56.187"	981.983

ANS: 1000 m in GRASS is 981.983 m according to proj4. May be due to datum change and not projection error? I'm not yet certain on this point.