Question about projection
Closed this issue · 9 comments
Hi Hesselbarth,
I have a question about which projection should I use to run the function calculate_lsm().
Using a wgs84 raster (2 classes, 0 and 1) I got the warming below, but the function worked and I got the result table.
"Warning message:
Caution: Coordinate reference system not metric - Units of results based on cellsizes and/or distances may be incorrect."
Then I reprojected the same raster to a system in meters (and I tried many possibilities: Albers Equal Area Conic, Equidistant Cylindrical, Lambert Azimuthal Equal Area, and the same projection of the object from package's memory "augusta_nlcd"...), but I got this error message for all projection systems:
"Erro: The area of the circumscribing circle is currently only implemented for equal resolutions."
What should I do in this case?
I am not the expert on CSR, but it seems like that the reprojection changes your resolution of cells differently in x and y direction.
What's the output of raster::res(x), where x is your landscape. Probably the values are not the same. Currently, we don't have a solution to calculate the circumscribing circle because it makes it really complicated to calculate its diameter.
So I think you need to find a metric CSR in which your cell resolution is identical in x and y direction. @Nowosad Any thoughts as CSR expert :)
Thank you for your answer!
I check it, but the resolution is the same...
[1] 0.00833 0.00833
Oh...maybe we need to relax our test then? Could you maybe send me your reprojected data either attached here or as an email to mhk.hesselbarth@gmail.com?
Hi Hesselbarth,
I just changed the command for albers projection and it worked now, thank you!
Now I used:
crs.albers <- CRS("+proj=aea +lat_0=-32 +lon_0=-60 +lat_1=-5 +lat_2=-42 +x_0=0 +y_0=0 +ellps=aust_SA +datum=WGS84 +units=m +no_defs")
The funny thing was that the result was the same in relation to the one that I performed before using WGS84.
Thank you very much
Hi @luaratourinho - map projections are a vast topic. Here I will mention just a few major things - hopefully, it could be useful for you.
- Geographic coordinate reference systems (CRS) have units in degrees. This is incompatible with how many landscape metrics are calculated. Algorithms behind many landscape metrics expect that one unit in one direction equals one unit in the other direction. However, this is not true for geographic CRS, where the length of one degree north is not only not equal to one degree east, but also one degree east on the equator has a different length than one degree on the mid-latitudes. See intro of https://geocompr.robinlovelace.net/reproj-geo-data.html.
- Therefore, if you want to calculate landscape metrics that depend on areas or distances, you need to reproject your data first. There are several hundreds of projections, and which one to use depends on the area you are studying, its size, etc.
- The
projectRaster()function of the raster package requires some customization to work with categorical rasters. See my examples below. Firstly, by default, it uses the bilinear interpolation, which changes the values of the original raster (it should only be used for continuous rasters). Secondly, I would encourage you to set some resolution to your output raster (which should be based on your knowledge about the data and your problem) and not to let the raster package decide on its resolution.
options(warn=-1)
library(raster)
#> Loading required package: sp
library(landscapemetrics)
# geographic CRS
podlasie_ccilc
#> class : RasterLayer
#> dimensions : 371, 457, 169547 (nrow, ncol, ncell)
#> resolution : 0.002777778, 0.002777778 (x, y)
#> extent : 22.23056, 23.5, 52.8, 53.83056 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
#> source : memory
#> names : ESACCI.LC.L4.LCCS.Map.300m.P1Y.1992_2015.v2.0.7
#> values : 10, 210 (min, max)
# resolution is in degrees
res(podlasie_ccilc)
#> [1] 0.002777778 0.002777778
hist(podlasie_ccilc)
# projected CRS #1 - wrong
podlasie_proj1 = projectRaster(podlasie_ccilc, crs = "EPSG:2180")
res(podlasie_proj1)
#> [1] 185 309
hist(podlasie_proj1)
# projected CRS #2 - discouraged
podlasie_proj2 = projectRaster(podlasie_ccilc, crs = "EPSG:2180",
method = "ngb")
res(podlasie_proj2)
#> [1] 185 309
hist(podlasie_proj2)
# projected CRS #2
podlasie_proj3 = projectRaster(podlasie_ccilc, res = c(300, 300),
crs = "EPSG:2180", method = "ngb")
res(podlasie_proj3)
#> [1] 300 300
hist(podlasie_proj3)
- Overall, I would also recommend you avoid using Proj4strings to set CRS. I explain why in this video - https://www.youtube.com/watch?v=Va0STgco7-4&list=PLf9p4wbX01At4CMbwCOx48mUPZVsF9Iw-&index=3. If possible - use EPSG codes instead.
- @mhesselbarth If I understand correctly, the original error
"The area of the circumscribing circle is currently only implemented for equal resolutions."means that the circumscribing circle algorithm expects a resolution in both dimensions (x and y) to be equal. Is this the expectation on many other metrics? If so - maybe it worth adding a check to thelandscape_check()function?
Sorry to post in a closed issue, but I noticed the warning message Coordinate reference system not metric - Units of results based on cellsizes and/or distances may be incorrect is no longer displayed when lsm_*() functions are called with landscapes in unprojected coordinates. The documentation for these functions also doesn't mention anywhere that landscapes should be projected, just saying results are in meters, which led me to believe that these functions were properly accounting for unprojected coordinates internally. Yet the areas and edge lengths returned by, for example, lsm_p_area() and lsm_c_te() for an unprojected landscape appear to be in units of degrees.
Do landscapemetrics functions "work" on unprojected coordinates and, if not, I'm wondering if the warning message should be re-instated and the documentation should state the results are in the units of the CRS rather than in meters as they currently state?
I will re-open the issue due to @mstrimas. That is something we need to look into since the warning message as well as some information in the docs should still be there. No idea where it went?
So technically all functions work for unprojected CRS, however, due to the units the ones based on areas and distances/lengths do not really make sense. For example, we use the cell resolution quite heavily (for example, the area is calculated by multiplying the number of cells times the resolution squared). We initially decided to only return a warning since many configuration metrics are not affected by this.
Oh this probably got lost quite a while ago when the internal structure was changed to use matrices whenever possible and not the raster itself. Thus, checking in each lsm_* function might be a bit tricky. Thus, we only check in calculate_lsm().
For now, I think adding some information to the documentation is the easiest fix.
I updated the documentation so could be closed again imho