0

I have a sf object with MULTIPOLYGONS that refers to flood episodes from 1800's onwards(https://data.gov.uk/dataset/16e32c53-35a6-4d54-a111-ca09031eaaaf/recorded-flood-outlines). Then I have another sf object, also MULTIPOLYGONS where I have information on local authorities for each year (1990-2019). What I want to do is to test - for each 1 year period - whether the each authority has been "touched" by the floods or not.

I have tried some things with sf's st_intersects, but I never quite got to take into consideration the time factor - what I usually get out is whether a local authority was ever affected by a flood.

dput of the flood object:

structure(list(name = c("06/11/2000_River Meon_Aerial(30)", "06/11/2000_River Meon_Aerial(32)", 
"Winter 2000/01gw_River Alre_Bishop's Sutton(2)", "Winter 2000/01gw_River Test_Braishfield(3)", 
"Winter 2000/01gw_River Test_Braishfield(2)", "Winter 2000/01gw_River Itchen_Bramdean(10)", 
"Winter 2000/01gw_River Itchen_Cheriton(3)", "07128FEG01369FEB01114", 
"06/11/2000_River Meon_Aerial(227)", "06/11/2000_River Meon_Aerial(209)"
), start_date = c("2000/11/06 00:00:00", "2000/11/06 00:00:00", 
"2000/12/01 00:00:00", "2000/11/01 00:00:00", "2000/11/01 00:00:00", 
"2000/12/01 00:00:00", "2000/12/01 00:00:00", "2050/01/01 00:00:00", 
"2000/11/06 00:00:00", "2000/11/06 00:00:00"), end_date = c("2000/11/07 00:00:00", 
"2000/11/07 00:00:00", "2001/03/01 00:00:00", "2000/12/01 00:00:00", 
"2000/12/01 00:00:00", "2001/03/01 00:00:00", "2001/02/01 00:00:00", 
"2050/01/01 00:00:00", "2000/11/07 00:00:00", "2000/11/07 00:00:00"
), Shape_Length = c(507.759112325055, 680.343081140914, 56.7885093275854, 
98.2313525071202, 28.7898295513574, 214.932542851196, 89.1577152757675, 
12.4839497962306, 26.2078119117354, 57.7726140088653), Shape_Area = c(5619.7706075058, 
19793.1513500502, 199.411550002538, 533.047249979913, 49.3870500007093, 
2015.71325006174, 505.290596007861, 12.2419819977947, 47.0836999942886, 
147.39440501066), geometry = structure(list(structure(list(list(
    structure(c(453881.42, 453890.2, 453911.45, 453902.89, 453901.3013, 
    453892.52, 453884.92, 453889.13, 453885.52, 453884.0264, 
    453882.22, 453879.07, 453876.99, 453876.23, 453875.87, 453877.3, 
    453877.3617, 453877.52, 453877.72, 453877.8, 453878.01, 453878.23, 
    453878.3, 453878.32, 453878.18, 453878.11, 453878.11, 453876.82, 
    453876.42, 453876.29, 453876.02, 453875.89, 453875.9, 453876.13, 
    453876.28, 453876.43, 453876.5, 453876.5068, 453877.1, 453877.98, 
    453879, 453880.02, 453881.11, 453881.47, 453883.14, 453865.09, 
    453855.13, 453862.75, 453866.2, 453860.78, 453856.3, 453841.3, 
    453835.92, 453836.65, 453845.21, 453858.06, 453867.02, 453865.68, 
    453872.55, 453881.42, 104299.59, 104288.57, 104296.57, 104280.810000001, 
    104277.8792, 104261.68, 104247.710000001, 104235.07, 104209.35, 
    104199.8727, 104188.41, 104165.689999999, 104149.880000001, 
    104141.199999999, 104135.52, 104133.42, 104133.3583, 104133.199999999, 
    104133.09, 104132.869999999, 104132.76, 104132.1, 104131.869999999, 
    104130.76, 104130.43, 104130.199999999, 104129.98, 104125.41, 
    104123.4, 104122.4, 104120.4, 104119.619999999, 104118.84, 
    104117.07, 104116.18, 104115.289999999, 104115.07, 104115.032500001, 
    104111.74, 104108.41, 104105.189999999, 104101.869999999, 
    104098.77, 104097.43, 104092.779999999, 104101.710000001, 
    104118.07, 104137.390000001, 104157.880000001, 104178.619999999, 
    104195.699999999, 104216.560000001, 104232.52, 104251.1, 
    104253.630000001, 104257.66, 104262.189999999, 104276.310000001, 
    104292.619999999, 104299.59), .Dim = c(60L, 2L)))), class = c("XY", 
"MULTIPOLYGON", "sfg")), structure(list(list(structure(c(454148.49, 
454168.39, 454169.4, 454169.42, 454168.24, 454167.1, 454067.28, 
454063.51, 454059.34, 454063.93, 454044.1, 453999.99, 453966.53, 
453988.76, 453992.4, 453996.77, 453995.32, 453993.56, 453999.97, 
454041.27, 454136.88, 454140.53, 454141.31, 454129.53, 454117.51, 
454118.71, 454117.28, 454131.32, 454128.48, 454148.49, 104302.15, 
104278.23, 104275.57, 104273.789999999, 104272.220000001, 104266.09, 
104284.289999999, 104280.92, 104282.1, 104294.27, 104298.720000001, 
104308.390000001, 104315.720000001, 104325.18, 104341.109999999, 
104360.4, 104377.84, 104398.73, 104418.48, 104421.35, 104430.140000001, 
104410.15, 104390.15, 104371.890000001, 104369.32, 104362.560000001, 
104357.09, 104347.67, 104330.51, 104302.15), .Dim = c(30L, 2L
)))), class = c("XY", "MULTIPOLYGON", "sfg")), structure(list(
    list(structure(c(460825.05, 460821.94, 460814.07, 460806.11, 
    460811.3, 460818.78, 460826.47, 460825.05, 131833.369999999, 
    131829.439999999, 131832.02, 131835.26, 131846.779999999, 
    131842.75, 131837.73, 131833.369999999), .Dim = c(8L, 2L)))), class = c("XY", 
"MULTIPOLYGON", "sfg")), structure(list(list(structure(c(437264.6, 
437261.43, 437257.92, 437253.16, 437249.16, 437247.33, 437245.63, 
437256.35, 437259.18, 437261.96, 437268.48, 437271.92, 437273.27, 
437274.35, 437272.76, 437264.6, 124097.550000001, 124089.970000001, 
124090.16, 124090.130000001, 124090.210000001, 124090.859999999, 
124092.41, 124114.609999999, 124120.529999999, 124122.77, 124124.039999999, 
124122.960000001, 124121.960000001, 124117.74, 124113.949999999, 
124097.550000001), .Dim = c(16L, 2L)))), class = c("XY", "MULTIPOLYGON", 
"sfg")), structure(list(list(structure(c(437362.11, 437359.25, 
437352.3, 437352.27, 437352.05, 437357.8, 437360.75, 437362.09, 
437362.11, 124349.449999999, 124347.1, 124347.720000001, 124351.6, 
124353.6, 124352.75, 124351.77, 124351.34, 124349.449999999), .Dim = c(9L, 
2L)))), class = c("XY", "MULTIPOLYGON", "sfg")), structure(list(
    list(structure(c(459329.38, 459296.31, 459261.92, 459254.01, 
    459261.94, 459281.8, 459305.57, 459326.74, 459347.94, 459353.19, 
    459329.38, 128039.119999999, 128031.18, 128035.119999999, 
    128040.369999999, 128052.35, 128056.25, 128061.529999999, 
    128061.550000001, 128058.9, 128053.619999999, 128039.119999999
    ), .Dim = c(11L, 2L)))), class = c("XY", "MULTIPOLYGON", 
"sfg")), structure(list(list(structure(c(458120.05, 458117.6, 
458116.06, 458112.54, 458111.28, 458107.7, 458104.4, 458101.54, 
458100.24, 458100.2436, 458100.78, 458104.31, 458111.8, 458120.41, 
458130.48, 458131.22, 458132.4315, 458132.96, 458133.29, 458132.42, 
458128.51, 458126.63, 458125.15, 458123.06, 458121.65, 458120.05, 
128587.85, 128587.710000001, 128587.689999999, 128587.76, 128587.75, 
128588.15, 128588.67, 128593.76, 128597.189999999, 128597.208799999, 
128599.98, 128603.02, 128604.220000001, 128605.32, 128606.99, 
128604.33, 128600.779100001, 128599.23, 128594.67, 128590.1, 
128589.17, 128588.810000001, 128588.68, 128588.33, 128588.09, 
128587.85), .Dim = c(26L, 2L)))), class = c("XY", "MULTIPOLYGON", 
"sfg")), structure(list(list(structure(c(458356.27, 458355.5, 
458354.73, 458354.09, 458353.65, 458353.5, 458353.65, 458354.09, 
458354.73, 458355.5, 458355.8421, 458356.27, 458356.91, 458357.0335, 
458357.35, 458357.5, 458357.35, 458356.91, 458356.27, 106768.75, 
106768.6, 106768.75, 106769.18, 106769.83, 106770.6, 106771.359999999, 
106772.01, 106772.439999999, 106772.6, 106772.528899999, 106772.439999999, 
106772.01, 106771.827500001, 106771.359999999, 106770.6, 106769.83, 
106769.18, 106768.75), .Dim = c(19L, 2L)))), class = c("XY", 
"MULTIPOLYGON", "sfg")), structure(list(list(structure(c(466662.36, 
466658.93, 466655.88, 466655.82, 466662.77, 466664.77, 466662.36, 
123515.810000001, 123515.32, 123517.949999999, 123522.07, 123522.039999999, 
123518.84, 123515.810000001), .Dim = c(7L, 2L)))), class = c("XY", 
"MULTIPOLYGON", "sfg")), structure(list(list(structure(c(465435.69, 
465435.2, 465433.59, 465431.97, 465429.79, 465425.6757, 465425.42, 
465425.0533, 465421.12, 465419.01, 465415.62, 465412.6, 465413.94, 
465421.58, 465428.86, 465435.69, 123960.890000001, 123954.76, 
123955.189999999, 123955.609999999, 123956.359999999, 123957.7817, 
123957.869999999, 123957.9791, 123959.15, 123959.57, 123960.970000001, 
123961.93, 123966.390000001, 123967.15, 123962.92, 123960.890000001
), .Dim = c(16L, 2L)))), class = c("XY", "MULTIPOLYGON", "sfg"
))), class = c("sfc_MULTIPOLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = 437245.63, 
ymin = 104092.779999999, xmax = 466664.77, ymax = 131846.779999999
), class = "bbox"), crs = structure(list(input = "OSGB 1936 / British National Grid", 
    wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n    BASEGEOGCRS[\"OSGB 1936\",\n        DATUM[\"OSGB 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n        BBOX[49.75,-9.2,61.14,2.88]],\n    ID[\"EPSG\",27700]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
10L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(name = NA_integer_, 
start_date = NA_integer_, end_date = NA_integer_, Shape_Length = NA_integer_, 
Shape_Area = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))

dput of the local authority: https://pastebin.com/4NveX0uU (had to limiti to 1st 5 rows due to size)

What I'm trying to do exactly is to check whether for every year (say 2005) each authority had or not had been affected from last May to May that year, so for 2005 I'm trying to see if any flood event intersected the authory boundaries from 15 May 2004 to 15 May 2005.

0

1 Answer 1

1

Please see example below. In the particular sample data, there seem to be no intersections with floods in any of the authorities in any year:

library(sf)

source("code_01_data.R")

# Data
dat1
## Simple feature collection with 10 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 437245.6 ymin: 104092.8 xmax: 466664.8 ymax: 131846.8
## projected CRS:  OSGB 1936 / British National Grid
##                                              name
## 1                06/11/2000_River Meon_Aerial(30)
## 2                06/11/2000_River Meon_Aerial(32)
## 3  Winter 2000/01gw_River Alre_Bishop's Sutton(2)
## 4      Winter 2000/01gw_River Test_Braishfield(3)
## 5      Winter 2000/01gw_River Test_Braishfield(2)
## 6      Winter 2000/01gw_River Itchen_Bramdean(10)
## 7       Winter 2000/01gw_River Itchen_Cheriton(3)
## 8                           07128FEG01369FEB01114
## 9               06/11/2000_River Meon_Aerial(227)
## 10              06/11/2000_River Meon_Aerial(209)
##             start_date            end_date Shape_Length
## 1  2000/11/06 00:00:00 2000/11/07 00:00:00    507.75911
## 2  2000/11/06 00:00:00 2000/11/07 00:00:00    680.34308
## 3  2000/12/01 00:00:00 2001/03/01 00:00:00     56.78851
## 4  2000/11/01 00:00:00 2000/12/01 00:00:00     98.23135
## 5  2000/11/01 00:00:00 2000/12/01 00:00:00     28.78983
## 6  2000/12/01 00:00:00 2001/03/01 00:00:00    214.93254
## 7  2000/12/01 00:00:00 2001/02/01 00:00:00     89.15772
## 8  2050/01/01 00:00:00 2050/01/01 00:00:00     12.48395
## 9  2000/11/06 00:00:00 2000/11/07 00:00:00     26.20781
## 10 2000/11/06 00:00:00 2000/11/07 00:00:00     57.77261
##     Shape_Area                       geometry
## 1   5619.77061 MULTIPOLYGON (((453881.4 10...
## 2  19793.15135 MULTIPOLYGON (((454148.5 10...
## 3    199.41155 MULTIPOLYGON (((460825 1318...
## 4    533.04725 MULTIPOLYGON (((437264.6 12...
## 5     49.38705 MULTIPOLYGON (((437362.1 12...
## 6   2015.71325 MULTIPOLYGON (((459329.4 12...
## 7    505.29060 MULTIPOLYGON (((458120 1285...
## 8     12.24198 MULTIPOLYGON (((458356.3 10...
## 9     47.08370 MULTIPOLYGON (((466662.4 12...
## 10   147.39441 MULTIPOLYGON (((465435.7 12...
dat2
## Simple feature collection with 5 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 418865.2 ymin: 506326.1 xmax: 478446.1 ymax: 537152
## projected CRS:  OSGB 1936 / British National Grid
##   objectid   lad15cd                       geometry year
## 1        1 E06000001 MULTIPOLYGON (((447213.9 53... 2016
## 2        2 E06000002 MULTIPOLYGON (((448489.9 52... 1994
## 3        3 E06000003 MULTIPOLYGON (((455834.1 52... 2001
## 4        4 E06000004 MULTIPOLYGON (((444157 5279... 2015
## 5        5 E06000005 MULTIPOLYGON (((423496.6 52... 2006
# Add 'year' column
dat1$year = substr(dat1$end_date, 1, 4)
dat1$year = as.numeric(dat1$year)

# Find intersections
for(i in unique(dat2$year)) {
    dat1_sub = st_union(dat1[dat1$year == i, ]) #to test if any district was affect by any flood that year
    dat2_sub = dat2[dat2$year == i, ]
    int = st_intersects(
        dat2_sub, 
        dat1_sub,
        sparse = FALSE
    )
    if(any(int)) dat2$flood[dat2$year == i] = int[ ,1] else dat2$flood[dat2$year == i] = FALSE
}

# Result
dat2
## Simple feature collection with 5 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 418865.2 ymin: 506326.1 xmax: 478446.1 ymax: 537152
## projected CRS:  OSGB 1936 / British National Grid
##   objectid   lad15cd                       geometry year
## 1        1 E06000001 MULTIPOLYGON (((447213.9 53... 2016
## 2        2 E06000002 MULTIPOLYGON (((448489.9 52... 1994
## 3        3 E06000003 MULTIPOLYGON (((455834.1 52... 2001
## 4        4 E06000004 MULTIPOLYGON (((444157 5279... 2015
## 5        5 E06000005 MULTIPOLYGON (((423496.6 52... 2006
##   flood
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.