How to find scale error for map projection?

Question

I've just upgraded from proj version 7.1 to 8.2, and some code has stopped working. I used to use proj to determine the map projection scale error at a point with the following code (source: https://lists.osgeo.org/pipermail/proj/2017-November/007496.html):

echo "-55.79072944 59.20415790" | proj -VS +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

Now when I run that, I get this error:

Rel. 8.2.0, November 1st, 2021
<proj>: 
can't initialize operations that take non-angular input coordinates
program abnormally terminated

This error is discussed on the proj GitHub: https://github.com/OSGeo/PROJ/issues/1619, but I don't see a solution from the discussion there - they suggest cct, but cct does not provide an equivalent for the -S flag: Causes estimation of meridional and parallel scale factors, area scale factor and angular distortion, and maximum and minimum scale factors to be listed.

I have a CSV file with 1000s of points and I'd like to estimate area error for each of them.

Is there a replacement for...

cat latlon.csv | proj -VS ${projstr} | grep Areal ...

?

Answer 1

In GRASS GIS with the latest proj, the following code works, and shows that the EPSG:3413 projection has a -6 % area error in N. Greenland and a +8 area error in S. Greenland.

GRASS> g.proj -j | grep -v crs
+proj=stere
+lat_0=90
+lat_ts=70
+lon_0=-45
+x_0=0
+y_0=0
+datum=WGS84
+units=m
+no_defs


GRASS> PROJSTR=$(g.proj -j | grep -v crs)

GRASS> echo "40 82.5" | proj -VS ${PROJSTR}
#Stereographic
#   Azi, Sph&Ell
#   lat_ts=
# +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84
# +units=m +no_defs +break_cs2cs_recursion +ellps=WGS84 +towgs84=0,0,0
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257224
#  squared eccentricity: 0.006694379990
Longitude: 40dE [ 40 ]
Latitude:  82d30'N [ 82.5 ]
Easting (x):   810474.21
Northing (y):  -70907.31
Meridian scale (h) : 0.97402441  ( -2.598 % error )
Parallel scale (k) : 0.97402441  ( -2.598 % error )
Areal scale (s):     0.94872355  ( -5.128 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 85d [ 85.00000000 ]
Max-min (Tissot axis a-b) scale error: 0.97402 0.97402



GRASS> echo "45 60" | proj -VS ${PROJSTR}
#Stereographic
#   Azi, Sph&Ell
#   lat_ts=
# +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84
# +units=m +no_defs +break_cs2cs_recursion +ellps=WGS84 +towgs84=0,0,0
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257224
#  squared eccentricity: 0.006694379990
Longitude: 45dE [ 45 ]
Latitude:  60dN [ 60 ]
Easting (x):   3323160.27
Northing (y):  -0.00
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 : 90d [ 90.00000000 ]
Max-min (Tissot axis a-b) scale error: 1.03943 1.03943