Dear GRASS List,
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the scale error for a given cell?
My understanding is that for a given projection, the "lat_ts" term of the proj4 string describing the projection is where 1 m in GRASS is actually 1 m. Elsewhere, errors grow. I'd like to track these.
-k.
Ken Mankoff wrote
Dear GRASS List,
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?My understanding is that for a given projection, the "lat_ts" term of the
proj4 string describing the projection is where 1 m in GRASS is actually 1
m. Elsewhere, errors grow. I'd like to track these.-k.
_______________________________________________
grass-user mailing list
grass-user@.osgeo
Doesn't it depends on the projection you're using?
Looking at e.g.
there is no lat_ts factor as this is an area preserving projection.
AFAIR there is now an area function in r.mapcalc in GRASS trunk to calculate
the pixel area.
HTH
-----
best regards
Helmut
--
Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Users-f3884509.html
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide
the scale error for a given cell?Doesn't it depends on the projection you're using?
Yes it does.
AFAIR there is now an area function in r.mapcalc in GRASS trunk to
calculate the pixel area.
Yes! I'm now motivated to upgrade to grass 7.3 to access the area() function.
Thank you for letting me know about this,
-k.
Dear Helmut and GRASS List,
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?My understanding is that for a given projection, the "lat_ts" term of
the proj4 string describing the projection is where 1 m in GRASS is
actually 1 m. Elsewhere, errors grow. I'd like to track these.AFAIR there is now an area function in r.mapcalc in GRASS trunk to
calculate the pixel area.
Is the scale error only in longitude and never in latitude? Or can it be in both? Is there a way to find the error for a length error rather than an area error? Specifically, for a vector line with an arbitrary heading?
Thanks,
-k.
Dear Ken,
On Thu, Nov 9, 2017 at 12:30 PM, Ken Mankoff <mankoff@gmail.com> wrote:
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?My understanding is that for a given projection, the "lat_ts" term of
the proj4 string describing the projection is where 1 m in GRASS is
actually 1 m. Elsewhere, errors grow. I'd like to track these.AFAIR there is now an area function in r.mapcalc in GRASS trunk to
calculate the pixel area.Is the scale error only in longitude and never in latitude? Or can it be in both?
Is there a way to find the error for a length error rather than an area error?
Specifically, for a vector line with an arbitrary heading?
This is probably a question for the PROJ mailing list?
Markus
Dear Helmut,
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?AFAIR there is now an area function in r.mapcalc in GRASS trunk to calculate
the pixel area.
There is. I've just built GRASS 7.5. But the area() function just returns nsres*ewres for each cell, it does not consider projection errors.
I'll look into the lower-level proj code for this.
Thanks,
-k.
On 20/11/17 21:20, Ken Mankoff wrote:
Dear Helmut,
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?AFAIR there is now an area function in r.mapcalc in GRASS trunk to calculate
the pixel area.There is. I've just built GRASS 7.5. But the area() function just returns nsres*ewres for each cell, it does not consider projection errors.
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. Comparing the output should give you a notion of projection error.
For visualization you might want to try the Indicatrix mapper Plugin in QGIS...
Moritz
On Tue, Nov 21, 2017 at 9:26 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:
On 20/11/17 21:20, Ken Mankoff wrote:
Dear Helmut,
On 2017-11-02 at 07:17, Helmut Kudrnovsky <hellik@web.de> wrote:
Ken Mankoff wrote
Does GRASS (or GDAL?, proj4?) have a function (or addon) to provide the
scale error for a given cell?AFAIR there is now an area function in r.mapcalc in GRASS trunk to calculate
the pixel area.There is. I’ve just built GRASS 7.5. But the area() function just returns nsres*ewres for each cell, it does not consider projection errors.
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.
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.
Markus M
Comparing the output should give you a notion of projection error.
For visualization you might want to try the Indicatrix mapper Plugin in QGIS…
Moritz
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
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@gmail.com> wrote:
On Tue, Nov 21, 2017 at 9:26 AM, Moritz Lennert <
mlennert@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.