Create routine to process TEOW Ecoregions data

[Jo-Schie]

For Teow Ecoregions we need a processing routine that intersects Ecoregions layer (polygon) with the PA polygons from WDPA. We need to calculate:

• Type of ecoregions that intersects with the PA
• Total area of each ecoregion in sqkm that intersects with the PA
• Share of each ecoregion area as part of total PA area.

The data can be assessed via API, see here.

Most probably for this issue we need to think a bit on how to automatically find an adequate equal-area projection system for the processed PAs? Any suggestions? Maybe @goergen95 has also an idea on this since this variable is also part of the old PA database, right @goergen95?

[Jo-Schie]

please also create a processing overview like in the pptx I send you.

[goergen95]

Hi!
So for the Latin America application, I think we calculated the UTM zone for each polygon based on its centroid.
For some edge cases, e.g. polygons with very large East-West extent, this actually would lead to bad results.
I would agree with @Jo-Schie that using some Lambert Azimuthal Equal Area projection based on the centroid of a polygon would deliver better results than blindly applying UTM zones without further checking the polygon. Is there already a script touching on this issue?

[Jo-Schie]

Hi @goergen95 and thanks for your response. I created such a script some time ago in the Accessibility Repo. However nobody ever reviewed it. Wanna take a look?

https://github.com/mapme-initiative/mapme.accessibility/blob/master/R/acc_areaproj.R

[goergen95]

Looks good to me! I am not sure if it is always a good choice for rasters because reprojecting will also change cell values. But I don't see an issue for projecting vector geometries this way when you are interested in area information.
Also, would you mind to explain how you set up the logic for the false easting/northing? How did you come up with the number 111000 in this line?

[goergen95]

Also, if I am not mistaken parameters +a and +b refer to the reference ellepsoid. Given the manual page of PROJ on LAEA I think it would be easier to just use the +ellps parameter and maybe also allow the user to change its default value?

[goergen95]

And one last thing: I think the function will fail when my_input is an sf-object. Is this behavior intended?

[Jo-Schie]

Thanks @goergen95 for looking at it. Before I adress your questions I have another doubt to solve. You send me a link a couple of days ago with standards for developing spatial software in R. I just went through that link again and it states the following:

The following standard applies to software which directly or indirectly relies on geographic data which uses or relies upon coordinate reference systems.

SP2.4 Geographical Spatial Software should be compliant with version 6 (and, ideally 7) of PROJ, and with wkt2 representations. The primary implication, described in detail in the articles linked to above, is that:

SP2.4a Software should not accept so-called “PROJ4-strings” previously used to specify coordinate reference systems.

Do I understand it correctly that prof4strings should not be used anymore? Because my function does just that...creating a proj4string. If that is the case, then we can just skip this discussion and I will have to think of an alternative anyways. I'm not sure, however, what that would really be.

[goergen95]

Yes, you are absolutely right. I fail at understanding how to easily customize CRSs with PROJ >4 so I tend to ignore the fact. Sorry for that. I will link a few resources that hopefully can guide us to a sustainable solution since proj4strings will no longer be supported. This manual explains how PROJ8 handles coordinate operations in general. A more R specific standpoint is available in a vignette of rgdal. As you might know, these changes to the PROJ API are quite recent (the r-spatial community announced to follow changes in GDAL and PROJ less than a year ago) and somewhat radical in breaking with the previous handling of transformations.

My take-away from these resources is that by dynamically searching for optimal conversions between CRSs and not relying on WGS84 as the central conversion hub in between, better accuracies are obtained. However, to achieve this PROJ relies on a database of authority compliant projection definitions. As I see it we either have the choice to select one of the pre-defined CRSs or we will have to customize the WKT definition to our own needs. Until now I have not been able to find a good resource explaining how to customize WKT definitions, but I guess I am missing something there. In fact, I have not been able to find a resource ever using a WKT definition that was not initialized by some EPSG code, which makes me wonder if maybe working with custom projections is also discouraged by recent changes in PROJ, at least implicitly.

[goergen95]

Here's a SO conversation also touching on custom LAEA projections with WKT.

[Jo-Schie]

Thanks, Darius, for these great resources. I just researched a bit on this topic as well with your ressources and could not find an easy solution. I think we have to develop this on the go and if we cannot find a solution maybe ask Oms supervisor at some point in time.

[Jo-Schie]

just for sake of completeness. Thanks for the link Darius that says that for authority free crs defintions proj4 is still valid. So I guess my function can be used.

To respond to your questions then:

1. How did you come up with the number 111000 in this line?
2. I think it would be easier to just use the +ellps parameter and maybe also allow the user to change its default value?
3. I think the function will fail when my_input is an sf-object. Is this behavior intended?

as to 1): I'm unsure why this is there. I copied the routine from here.
as to 2): Also unsure. The whole idea of this is to have an automated definition for the CRS so in Theory there should be no user but the CRS should be defined on the extent of another spatial object. So I would guess it does not make sense to provide that option, at least without an automated routine to also generate the best ellipsoid...
as to 3): Correct and bad practice after consulting the link you send some time ago. I will have to update the accessibility package and I put this in the repo as a new issue.

[goergen95]

FYI, the rOpenSci has asked r-spatial to review their standards for spatial software. So these are likely to be changed a little bit in the near future. Also, note this comment on SP2.4a stating that usage of proj4 strings is still valid for some conversions.