[PROJ] Correcting map projection errors

[daan]

It’s simple enough to just implement yourself. Use the authalic latitude to project your geodetic latitude to the sphere, and then use the spherical cylindric equal-area to get to the plane. From there you’re dealing with rectangles.

Authalic latitude: 
  mathworld.wolfram.com/AuthalicLatitude.html

Lambert: x= longitude, y = sin(latitude)

Cheers,
— daan

> On Apr 23, 2019, at 10:55, Ken Mankoff <mankoff at 3m411.com> wrote:
> 
> 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.
>