<div dir="ltr"><div dir="ltr">daan wrote:<br>> Just project your cells using ellipsoidal cylindric equal area<br>> projection and compute the area of the resulting rectangles. This<br>> gives you perfect accuracy at the minimum possible computational<br>> effort.<br><br>Ken replies:<br>> This sounds like it might work, but I'm not sure how to do this,<br>> i.e. what the correct EPSG code or proj string is for this<br>> projection.<br><br>While this may work, and provide an acceptable accuracy, it is certainly not recommendable.<br><br>Essentially, it translates into resampling into a different coordinate system, typically by linear interpolation.<br><br>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!<br><br>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).<br><br>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.<br><br>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.<br><br>Ken says: <br>> I have various products that cover all of Greenland at<br>> 200 x 200 m. This is 100 million cells. I'd like to estimate<br>> the error for each cell, 1x. For each individual one I would<br>> like to say "this is really 40100 m^2" and<br>> "this is really 39234 m^2".<br><br>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.<br><br>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.<br><br>/Thomas<br><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Den tir. 23. apr. 2019 kl. 20.03 skrev daan <<a href="mailto:strebe@aol.com">strebe@aol.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">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.<br>
<br>
Authalic latitude: <br>
  <a href="http://mathworld.wolfram.com/AuthalicLatitude.html" rel="noreferrer" target="_blank">mathworld.wolfram.com/AuthalicLatitude.html</a><br>
<br>
Lambert: x= longitude, y = sin(latitude)<br>
<br>
Cheers,<br>
— daan<br>
<br>
> On Apr 23, 2019, at 10:55, Ken Mankoff <<a href="mailto:mankoff@3m411.com" target="_blank">mankoff@3m411.com</a>> wrote:<br>
> <br>
> Hi Thomas and Daan,<br>
> <br>
> Daan wrote:<br>
>> Just project your cells using ellipsoidal cylindric equal area<br>
>> projection and compute the area of the resulting rectangles. This<br>
>> gives you perfect accuracy at the minimum possible computational<br>
>> effort.<br>
> <br>
> 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.<br>
> <br>
> <br>
> <br>
> On 2019-04-23 at 19:00 +0200, Thomas Knudsen <<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>> <br>
>> I think I understand your case better now: Your base data are gridded<br>
>> (either model output or EO), and you want to know the “real” area of a<br>
>> cluster of grid cells each nominally being 200 m x 200 m.<br>
> <br>
> Thomas - yes I'm sorry if I didn't explain my model use case clearly.<br>
> <br>
> 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".<br>
> <br>
> 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<br>
> <br>
> "echo ${lon} ${lat} | proj -VS ${PROJSTR}"<br>
> <br>
> 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)<br>
> <br>
> 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.<br>
> <br>
> 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...<br>
> <br>
> 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.<br>
> <br>
>  -k.<br>
> <br>
<br>
</blockquote></div>