Search code examples
rasterr-rasterr-spr-sf

Zonal operation using an overlapped-ratified raster which comes from two other ratified rasters

I have the following reclassify & ratified rasters that I am trying to intercept/overlay/overlap to get a new one. The idea is to get a new overlayed raster from the interception areas of the rasters R1 and R2. Once done this, I would do zonal operations. Here the R1, R2, ED rasters.

R1:
class      : RasterLayer 
dimensions : 1399, 1855, 2595145  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : UK_GDP_2010_PPP_percapita_km2 
values     : 1, 6  (min, max)
attributes :
       ID AG
 from:  1  a
  to :  6  f

R2:
class      : RasterLayer 
dimensions : 1399, 1855, 2595145  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 1, 5  (min, max)
attributes :
 ID AG
  1  A
  2  B
  3  C
  4  D
  5  E

Here the code to intercept/overlay/overlap

1st approach
library(sf)
library(tidyverse)
R1_SPDF <- as(R1,'SpatialPolygonsDataFrame')
R1_SPDF <- st_as_sf(R1_SPDF)

R2_SPDF <- as(R2,'SpatialPolygonsDataFrame')
R2_SPDF <- st_as_sf(R2_SPDF)

R3 <- st_intersection(R1_SPDF, R2_SPDF)

R3:
Simple feature collection with 1174501 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -8.166667 ymin: 50.01667 xmax: 1.541667 ymax: 59.55
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
    UK_GDP_2010_PPP_percapita_km2 layer
1                               1     1
2                               1     1
2.1                             1     1
3                               1     1
3.1                             1     1
4                               1     1
4.1                             1     1
1.1                             1     1
5                               1     1
6                               1     1
                      geometry
1   POLYGON ((-1.641667 59.55, ...
2   POLYGON ((-1.633333 59.55, ...
2.1 POLYGON ((-1.633333 59.55, ...
3   POLYGON ((-1.625 59.55, -1....
3.1 POLYGON ((-1.625 59.55, -1....
4   POLYGON ((-1.616667 59.55, ...
4.1 POLYGON ((-1.616667 59.55, ...
1.1 POLYGON ((-1.641667 59.5416...
5   POLYGON ((-1.65 59.54167, -...
6   POLYGON ((-1.641667 59.5416...

However, I am not sure if this result is what I am looking for, because I expect a new ratified raster R3 with the overlayed areas formed with the combination/interception of R1 and R2 areas (R3 areas ratify: Aa,.. Af,…,Ea,…Ef) or something like that. 

Expected R3:
class      : RasterLayer 
dimensions : #from intersection 
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : layer 
values     : 1, 30  (min, max) #aproximately 30 because R1: ID=6, and R2: ID=5.
attributes :
 ID AGnew
  1  Aa
  2  Ab
  .  .
  .  .
  .  .
  30 Ef

Here a second try using the raster package:

2nd approach
library(raster)
R3.1 <- intersect(R1_SPDF, R2_SPDF)
R3.1: 
Simple feature collection with 485611 features and 0 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -8.65 ymin: 49.875 xmax: 1.766667 ymax: 60.84167
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
                     geometry
1  POLYGON ((-0.9 60.84167, -0...
2  POLYGON ((-0.8916667 60.841...
3  POLYGON ((-0.8833333 60.841...
4  POLYGON ((-0.875 60.84167, ...
5  POLYGON ((-0.85 60.84167, -...
6  POLYGON ((-0.8416667 60.841...
7  POLYGON ((-0.9 60.83333, -0...
8  POLYGON ((-0.8916667 60.833...
9  POLYGON ((-0.8833333 60.833...
10 POLYGON ((-0.875 60.83333, ...

Once I got the R3, I expect to do the following zonal operation. Sum the values of ED raster within the R3 overlayed areas. 

sum_R <- zonal(ED, R3, "sum")

Any recommendation is very welcome.

Solution

  • I think the answer in not clear at all, though I also think you already have the answer. Perhaps your main issue here is that you are working with relative large datasets and you do not manage to see what is going on at every stage.

    So for the sake of simplicity I propose a smaller instance of your problem. Notice that when asking in stackoverflow this approach may be more useful to get help. As a matter of fact your files may be too heavey to download and use in some computers, avoiding people to get involved.

    But let's go to the code. 

    library(tidyverse)
library(raster)
library(sf)
library(tmap)

    You basically have two rasters with the same resolution, location and extent as the followings: 

    R1 <- raster(ncol=10, nrow=10, xmn = 0, xmx=10, ymn = 0, ymx = 10)
R2 <- raster(ncol = 10, nrow = 10, xmn = 0, xmx=10, ymn = 0, ymx = 10)
values(R1) <- c(rep(1,ncell(R1)/2), rep(2,ncell(R1)/2))
values(R2) <- rep(10,ncell(R2))

    with these kind of files you can perfomer directly many operations. The same extent and resolution is a plus: 

    Rsum <- R1 + R2
plot(R1) #See the legend of the plot
plot(R2) #See the legend of the plot
plot(Rsum) #See the legend of the plot

    I do not see why you need to perform many other operation to get to a zonal operation.

    It would make sense if your rasters were different. For example: 

    R2_Alt <- raster(ncol = 5, nrow = 2, xmn = 0, xmx=10, ymn = 0, ymx = 10) #Alt for alternative
values(R2_Alt) <- rep(10,ncell(R2))

    In the case of using R2_Al a simple operation could produce a mistake:

    Test <- R2_Alt + R1

    So it make sense to move to vector files: 

    R1_sp <- as(R1, "SpatialPolygonsDataFrame")
R2_Alt_sp <- as(R2_Alt, "SpatialPolygonsDataFrame")

R1_sf <- st_as_sf(R1_sp)
R2_Alt_sf <- st_as_sf(R2_Alt_sp)

names(R1_sf)[1] <- "V1" #Just to get differente variable names
names(R2_Alt_sf)[1] <- "V2_A"


tm_shape(R1_sf) + tm_polygons(border.col = "black") +
 tm_shape(R2_Alt_sf) + tm_polygons(border.col = "red", lwd = 4, alpha = 0)

    enter image description here

    And here is where having a smaller instance of your problem could shed light on your issues: 

    Int <- st_intersection(R2_Alt_sf, R1_sf)
View(Int)
plot(st_geometry(Int))

    You can also see the outcome with the former R2 raster:

    R2_sp <- as(R2, "SpatialPolygonsDataFrame")
R2_sf <- st_as_sf(R2_sp)
names(R2_Alt_sf)[1] <- "V2_A"

Int2 <- st_intersection(R2_sf,R1_sf)

    Lastly, be sure that the operation you need is intersection. see This.

    I know this answer is not really an answer, but I hope that it helps you to get closer. Let me know what you think and put here all the feedback.

    all the best