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.

Sorry if this is too abstract, but I struggled with ways to paste or represent the data here. If anyone has suggestions for that, I'm happy to hear.

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.

Thanks!

4
  • The most straightforward approach is to use a for loop, going over each unique year, (1) taking a subset of both layers with the polygons for that year, (2) calculating a new column in the authorities (subset) layer indicating whether each authority was flooded, based on st_intersects and the (subset) of flood episode layer, (3) combining back to a complete authorities layer for all years including the new columns. If you can prepare a small reproducible example of both datasets it will be easier to answer with a specific code example. Dec 1 '20 at 9:44
  • Oh I see! Great. thanks a lot. I added 2 dput samples and added some more context of what I'm trying to do at the bottom. Hope it helps!
    – AntVal
    Dec 1 '20 at 10:54
  • Thanks! The second object (pastebin.com/NjxFGDsy) isn't available Dec 1 '20 at 13:58
  • Sorry. should work now: pastebin.com/4NveX0uU
    – AntVal
    Dec 1 '20 at 14:46
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
6
  • 1
    This did it! Thank you so much for all the help & patience with me!
    – AntVal
    Dec 2 '20 at 13:54
  • Great, happy to hear! Dec 2 '20 at 14:22
  • sorry to annoy you again, but there's a mistake. I wished to test the intersections for each lad15nmyear in dat1, and not to get them to have the same result for each year. As it stands right now, if any of the lad15nm values has int == 1, then ALL get it for that value of year. What I want is to have it specifically for each lad15nmyear pair. Sotty for not being clearer again!
    – AntVal
    Dec 10 '20 at 8:59
  • If I understand correctly, you want to find out for each authority whether it intersects with each flood, in each year. If so, the information cannot be contained in one layer as above, because the floods are different each year. In what way do you need the data to be arranged? For example, in a table with authority / flood / year / flooded(yes/no) columns? Dec 11 '20 at 12:14
  • That's it, I probably was not clear: I just want to check whether each authority i in year j is flooded by any of floods that occured in year j, but only those of year j. That's what I think I get from the small edit I did in your reply. Ideally I need only authority/flood/year/flooded (Y/N).
    – AntVal
    Dec 11 '20 at 12:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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