[Proj] Pixel area errors

Mon Nov 20 12:56:56 PST 2017

Ken,

lat_ts is a setup parameter used in a few of the projections supported by
proj (cea, merc, stere, wag3 and wink1) - it is not a magic incantation you
can slap on any projection definition string to force it to be "exact"
along any given latitude.

But the proj program has functionality to provide a large amount of
projection distortion information: Try running

echo 12 55  |  proj -VS +proj=merc +ellps=GRS80 +lat_ts=55

Which will give you information about the distortion of a point at 12 E, 55
N (near Copenhagen, Denmark), for a mercator projection with a latitude of
true scale of 55. Not surprising, the scale is 1:

#Mercator
# Cyl, Sph&Ell
# lat_ts=
# +proj=merc +ellps=GRS80 +lat_ts=55
#--- following specified but NOT used
# +ellps=WGS84
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 12dE [ 12 ]
Latitude:  55dN [ 55 ]
Easting (x):   767929.55
Northing (y):  4211972.20
Meridian scale (h) : 1.00000000  ( 8.866e-009 % error )
Parallel scale (k) : 1.00000000  ( -1.066e-010 % error )
Areal scale (s):     1.00000000  ( 8.759e-009 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ -0.00000000 ]
Max-min (Tissot axis a-b) scale error: 1.00000 1.00000

Then try changing +lat_ts=55 to +lat_ts=54:

#Mercator
# Cyl, Sph&Ell
# lat_ts=
# +proj=merc +ellps=GRS80 +lat_ts=54
#--- following specified but NOT used
# +ellps=WGS84
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 12dE [ 12 ]
Latitude:  55dN [ 55 ]
Easting (x):   786909.29
Northing (y):  4316073.03
Meridian scale (h) : 1.02471546  ( 2.472 % error )
Parallel scale (k) : 1.02471546  ( 2.472 % error )
Areal scale (s):     1.05004178  ( 5.004 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ -0.00000000 ]
Max-min (Tissot axis a-b) scale error: 1.02472 1.02472

Indeed a quite useful feature of proj! That said, you will probably also
need some basic literature about projections to make real use of this
feature.

/thomas

2017-11-20 20:28 GMT+01:00 Micah Cochran <mcochran at athensal.us>:

> Hello Ken,
>
> Just for completeness a projection is a flat representation of earth.  The
> earth is usually represented as an oblate spheroid.  Any projection has
> some level of approximation between distance, area, direction, and form
> (does the landmass look like what it should look like).  Some projections
> are good used for surveying application, others are better for representing
> a really large area.
>
> For basic information about projections:
> http://axismaps.github.io/thematic-cartography/articles/projections.html
>
> If you limit the scope of area, you can very accurate measurements.
>
> It depends upon what projection that you are using to what kinds of errors
> you will get.  Mercator projection is very accurate at the 0 deg latitude.
> Transverse Mercator project is very accurate at along a line of longitude.
>
> Errors can be both latitude and longitude, depending on the projection.
> For example, a polar projection will have distortion anywhere but the
> center (the pole).
>
> Tissot's indicatrix is used as a graphic illustration of a projection's
> distortion.
> https://en.wikipedia.org/wiki/Tissot%27s_indicatrix
>
> You can calculate a grid distance, which can be very close or not close at
> all based on the projection.    Geodetic/geodesic distance algorithms from
> the GeographicLib::Geodesic libary (included with PROJ.4) has very accurate
> distance calculations.
> http://proj4.org/geodesic.html
>
> I'm not to sure about keeping track of the errors of a particular
> projection or how to exactly answer your specific questions.  I hope that
> gives you some pointers.
>
> In a slightly different tangent about using US Customary units for area
> and exactitude, I was looking at the difference between calculating acres
> using international feet (0.3048 m) and US Survey feet (0.3048006 m).  The
> size of this City is about 1000 sq km.  When converting to acres (area)
> using the two different feet as the basis of acres, I got a difference of
> 100 acres (40 hectares or 0.4 sq km) between the asnwers.  US Customary
> units make me crazy.
>
> Micah Cochran
>
> On Mon, Nov 20, 2017 at 11:13 AM, Ken Mankoff <mankoff at gmail.com> wrote:
>
>> Dear List,
>>
>> I've just learned that GRASS has a function to provide the true scale for
>> a given cell. From this, I can determine the map projection error.
>>
>> 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.
>>
>> I asked the following on the GRASS list and was advised that this list
>> might be a better place to ask:
>>
>> 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?
>>
>> Thank you,
>>
>>   -k.
>> _______________________________________________
>> Proj mailing list
>> Proj at lists.maptools.org
>> http://lists.maptools.org/mailman/listinfo/proj
>>
>
>
>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20171120/ffc3683c/attachment.html>