CODEWITHSUNDEEP
rspatial
jimken
jimken

Reputation: 117

Calulate area of cells in lat lon grid in r

I want to create a spatial grid 2x1 degrees and calculate the area in each grid cell. I used the following code but the answer is clearly not correct with cell size calculations coming out at 4500-14 000 km2. Can anyone point out where I am going wrong.

library(sf)
library(raster) 
library(dplyr)

xlim <- c(-45, 27)
ylim <- c(53, 80)

#Create grid
grid <- data.frame(lat = ylim, lon = xlim) %>% 
  st_as_sf(coords=c('lon','lat'), crs=4326) %>% 
    st_make_grid(cellsize = c(2,1)) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(cellid = row_number()) 

grid$area <- st_area(grid) #area of each "square"
grid$area_km2 <- units::set_units(grid$area, "km^2") # in km2

Upvotes: 0

Views: 525

Answers (1)

Chris
Chris

Reputation: 2286

Your cells are going to have different areas. This is easiest to demonstrate by changing ylim = c(-40, 40) equal above and below equator, depending on where you're standing, run your code above, then, after grid$area_km2 <-:

grid_km2_rle <- rle(as.vector(grid$area_m2))
head(grid4_km2_rle)
$lengths
 [1] 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36
[26] 36 36 36 36 36 36 36 36 36 36 36 36 36 36 72 36 36 36 36 36 36 36 36 36 36
[51] 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36
[76] 36 36 36 36

$values
 [1] 19098.63 19366.15 19627.69 19883.17 20132.53 20375.71 20612.63 20843.24
 [9] 21067.48 21285.29 21496.60 21701.38 21899.56 22091.10 22275.95 22454.06
[17] 22625.38 22789.88 22947.52 23098.25 23242.04 23378.85 23508.65 23631.41
[25] 23747.11 23855.70 23957.17 24051.49 24138.65 24218.61 24291.36 24356.88
[33] 24415.17 24466.19 24509.95 24546.44 24575.63 24597.54 24612.14 24619.44
[41] 24612.14 24597.54 24575.63 24546.44 24509.95 24466.19 24415.17 24356.88
[49] 24291.36 24218.61 24138.65 24051.49 23957.17 23855.70 23747.11 23631.41
[57] 23508.65 23378.85 23242.04 23098.25 22947.52 22789.88 22625.38 22454.06
[65] 22275.95 22091.10 21899.56 21701.38 21496.60 21285.29 21067.48 20843.24
[73] 20612.63 20375.71 20132.53 19883.17 19627.69 19366.15 19098.63

so, -40 and 40 have the same area:

min(grid_km2_rle$values)
[1] 19098.63 
# whereas at equator
grid_km2_rle$values[which.max(grid_km2_rle$values)]
[1] 24619.44

effects of sphere to plane.

Upvotes: 1

Related Questions