Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Misleading area calculations using occurrence data (AOOarea function) #3

Open
gepinillab opened this issue Mar 3, 2023 · 0 comments

Comments

@gepinillab
Copy link
Member

Hi again,

When calculating the AOO based on the occurrence of a particular feature, it is important to consider that if multiple occurrences of that feature are present in the same pixel, the area calculation will be skewed.

To illustrate this point, I have created an example using occurrence data in Manhattan. Using a dataset of 1000 random points within Manhattan, I attempted to estimate the area of Manhattan using these points. However, because several points are in the same pixels, the resulting area estimation was over 4000 km², which is clearly an inaccurate and misleading result. Manhattan's actual area is only 60 km².

I have included the code for the example in Manhattan below:

library(tigris)
library(terra)
library(changeRangeR)
library(raster)

# Get Manhattan polygon
mh <- tigris::tracts(state = '36', county = '061')
mh <- tigris::erase_water(mh)
sf::sf_use_s2(FALSE)
mh <- sf::st_union(mh) %>% sf::st_buffer(dist = 0)
plot(mh)
## Area
terra::expanse(terra::vect(mh), unit = "km")
### 59.1 km2 (match Wikipedia)

## To WGS 84
mh <- terra::project(terra::vect(mh), y = "EPSG:4326") %>% 
  sf::st_as_sf()
# Create 1000 points inside Manhattan
pts <- terra::spatSample(terra::vect(mh), 1000, method = "regular") %>% 
  sf::st_as_sf()
# AOO based on points ###
## cRR
### Points need to be in WGS84 and it requires a template raster
r <- rast(mh, res = 1/120)
values(r) <- 1
rWGS84 <- terra::mask(r, mh)
### Run
crrPts <- changeRangeR::AOOarea(methods::as(r, "Raster"), 
                                sf::st_coordinates(pts))
print(crrPts$area) # Bug: 4072 km2!!!

Let me know if you have any questions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant