[PROJ] Correcting map projection errors

[Ken Mankoff]

Hi Thomas and Daan,

Daan wrote:
> Just project your cells using ellipsoidal cylindric equal area
> projection and compute the area of the resulting rectangles. This
> gives you perfect accuracy at the minimum possible computational
> effort.

This sounds like it might work, but I'm not sure how to do this, i.e. what the correct EPSG code or proj string is for this projection. After some searching I *think* "Lambert Cylindrical Equal Area" might be what you're referring to, which has EPSG:9835, but GRASS does not support that.



On 2019-04-23 at 19:00 +0200, Thomas Knudsen <knudsen.thomas at gmail.com> 
> 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.

Thomas - yes I'm sorry if I didn't explain my model use case clearly.

I have various products that cover all of Greenland at 200 x 200 m. This is 100 million cells. I'd like to estimate the error for each cell, 1x. For each individual one I would like to say "this is really 40100 m^2" and "this is really 39234 m^2".

Right now I do this as I showed in my original email. For every 100 M cell center (x,y) I get the (lon,lat) from GRASS, pass that to proj with

"echo ${lon} ${lat} | proj -VS ${PROJSTR}"

and then capture the output, and have the area error for that cell. The result of this step is a new raster, resolution 200 x 200 m, 100 M cells. Cell values range from ~8 to ~-6 (EPSG:3413 has up to 8 % error)

Based on your suggestions that planimeter is better than proj with this stuff, I will rewrite the code to pass the outline of each individual cell to the planimeter tool, and then capture and parse the output to get a matching raster: 200 x 200 m covering my entire domain, with a scale factor at each cell location that represents the projection error.

I am now concerned that the GRASS step of "convert each cell corner from (x,y) to WGS84 (lon,lat)" may introduce its own errors...

But, back to my original question... The proj-generated error raster varies *very slowly*. I take this to mean that your original comment about "valid for infinitesimal items only" may not be an issue. I don't think there is a significant error introduced by my cell-center-point method. However, from the lengthy thread you linked to and your comments about planimeter, proj may not be providing correct results.

  -k.