[PROJ] Correcting map projection errors

[Thomas Knudsen]

Ken,

I think I understand your case better now: Your base data are gridded
(either model output or EO), and you want to know the “real” area of a
cluster of grid cells each nominally being 200 m x 200 m.

Hence, nominally, your cells are 40000 m2, but really, they differ. So I
presume your expected accuracy is +/- 1 cell, i.e. +/-40000 m2.

Doing the computation using the geodesic planimeter, you would need to
trace the outline of each cell cluster (corner coordinate of each outer
cell), convert to geographical WGS84, then computing the area of the
resulting tiny-sided polygon using the planimeter tool.

To get a feeling of the precision you may obtain, try comparing the
geodesic area of a cell expressed as a polygon defined by its corner
coordinates with the properly scaled value for your original approach:
Charles Karney states that the accuracy of his algorithm is on the order of
0.1 m2/vertice, i.e. 0.4 m2 for a polygon tracing the outline of a grid
cell of nominally 40000 m2 - i.e. a relative accuracy of approximately 1e-5.

Since the original mapping to the projected grid may be dubious, it may
even be a (line-of-code wise) fine shortcut to compute individual geodesic
areas for each raster cell, then summing up everything (I suppose the grid
parameters are static, and you will only have to do this once).

/Thomas


Den tir. 23. apr. 2019 kl. 11.47 skrev Ken Mankoff <mankoff at 3m411.com>:

> Hi Thomas,
>
> Thanks for the reply.
>
> > But note that the figures you obtain are valid for infinitesimal items
> > only. For large lengths/areas you should do the computations with
> > geodesics on the ellipsoid
>
> Can you clarify what you mean by this? Or more specifically, am I doing it
> incorrectly? I'm exporting a very large grid but at 200 x 200 m resolution
> to a CSV file, then using proj to calculate the error every 200 m. Is that
> close enough to infinitesimal? Or define "large lengths" please. I'm then
> using the errors to correct linear features at the ~200 m cell resolution.
>
> > using e.g. Charles Karney’s Planimeter tool over at
> > https://geographiclib.sourceforge.io/cgi-bin/Planimeter (online
> > service, includes link to source etc.)
>
> I've downloaded and compiled the code. Planimeter takes input in lat/lon,
> UTM/UPS, or MGRS. I'm working in EPSG:3413. I can provide lon,lat at each
> cell rather than x,y but I'm not sure if this introduces new errors.
>
> > Also, a day where you have too much time on your hand, see the
> > discussion over at
> >
> http://osgeo-org.1560.x6.nabble.com/Re-Qgis-user-New-Features-in-Shape-Tools-3-2-0-td5378898.html
>
> A nice read - not that I understood all of it - and one reason I use and
> contribute to OSS/FS.
>
>   -k.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190423/6c7468e5/attachment-0001.html>