[PROJ] Correcting map projection errors

[Thomas Knudsen]

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.

Ken replies:
> 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.

While this may work, and provide an acceptable accuracy, it is certainly
not recommendable.

Essentially, it translates into resampling into a different coordinate
system, typically by linear interpolation.

If the grid comes from model output, and grid values are point values (not
cell mean values), that may make sense. But if the model output is cell
mean values, you would (in loose terms) need to convolve the input grid
with a kernel the shape of the output grid cell, to get proper values for
the output grid. Point- and block support are geostatistically different!

If the grid is “level 3 style” EO data (and from the EPSG:3413 description,
I suppose that would mean Radarsat), the N Greenland data density will
probably be so high that regridding would make very little difference when
aggregating over larger areas (since it’s already heavily aggregated).

This, however, assumes that you are interested in mean values, and not (a
time series of) anomalies. Also, the data density in S Greenland will be
smaller, so down there, the resampling may have a higher (but probably
still small) influence on the areal mean.

So I still think resampling using authalic latitudes is a dubious idea. But
only a domain expert (in casu Ken) can have an informed opinion on whether
the deviation is acceptable.

Ken says:
> 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".

Yes, this makes sense. But remember to be careful when computing aggregate
statistics of large numbers of individual values: You may be able to save a
few additional significant digits by storing the anomaly (i.e. the area
minus 40000 m2), rather than the area herself.

But we’re drifting far outside the subject domain of the PROJ mailing list,
so Ken: I will contact you directly - I see from your personal home page
that you work for GEUS, which means we are part of the same DK branch of
government, so I could probably drop by for further discussion in person,
if needed.

/Thomas


Den tir. 23. apr. 2019 kl. 20.03 skrev daan <strebe at aol.com>:

> 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.
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190424/aa0d9c6c/attachment.html>