0

I have a map of Mexico and a data.frame with 300 weather stations. I drew 10km buffers around each station using st_buffer. The buffers are stored in buffers and the map is stored in mex.

I plotted the map and the buffers. My goal is to assign the station number of a given buffer to all the polygons that are touched by that buffer.

My map and buffers look like this.

The map can be downloaded from this site. The CRS has to be changed from 5332 to 4326 first.

These are the first two buffers:

> dput(buffers[1:2, ])
structure(list(station_number = 1004:1005, station_alt = c(1925, 
1844), month = c(9L, 9L), Mean_min = c(11.6, 12.75), Mean_max = c(26.9333333333333, 
26.85), months_observed = c(5L, 5L), geometry = structure(list(
    structure(list(structure(c(-102.193501864824, -102.188739639421, 
    -102.184005627476, -102.179312845387, -102.174674195097, 
    -102.170102428304, -102.165610111091, -102.161209589104, 
    -102.156912953352, -102.15273200674, -102.14867823144, -102.144762757164, 
    -102.14099633046, -102.137389285095, -102.133951513622, -102.130692440198, 
    -102.127620994743, -102.124745588492, -102.122074091026, 
    -102.119613808832, -102.117371465453, -102.115353183273, 
    -102.113564467008, -102.112010188909, -102.110694575746, 
    -102.109621197594, -102.108792958435, -102.108212088624, 
    -102.107880139208, -102.107797978126, -102.107965788291, 
    -102.108383067551, -102.109048630528, -102.109960612321, 
    -102.111116474066, -102.112513010319, -102.114146358253, 
    -102.116012008624, -102.11810481849, -102.12041902562, -102.122948264575, 
    -102.125685584399, -102.128623467869, -102.13175385226, -102.135068151548, 
    -102.138557280014, -102.142211677156, -102.146021333863, 
    -102.149975819769, -102.154064311712, -102.158275623224, 
    -102.162598234981, -102.167020326114, -102.171529806314, 
    -102.176114348636, -102.180761422931, -102.18545832979, -102.190192234944, 
    -102.194950204009, -102.199719237488, -102.204486305951, 
    -102.209238385282, -102.21396249192, -102.218645717996, -102.223275266261, 
    -102.227838484744, -102.232322901016, -102.236716255998, 
    -102.241006537204, -102.245182011349, -102.24923125622, -102.253143191741, 
    -102.256907110134, -102.260512705114, -102.263950100018, 
    -102.267209874814, -102.270283091899, -102.273161320625, 
    -102.275836660488, -102.278301762903, -102.280549851515, 
    -102.282574740994, -102.284370854233, -102.285933237938, 
    -102.287257576532, -102.288340204344, -102.28917811605, -102.289768975321, 
    -102.290111121663, -102.290203575419, -102.290046040907, 
    -102.289638907696, -102.288983250005, -102.288080824209, 
    -102.286934064475, -102.285546076523, -102.283920629526, 
    -102.28206214616, -102.279975690846, -102.277666956193, -102.275142247682, 
    -102.272408466631, -102.269473091493, -102.266344157511, 
    -102.263030234815, -102.259540405001, -102.255884236252, 
    -102.25207175709, -102.248113428804, -102.244020116654, -102.239803059918, 
    -102.235473840866, -102.231044352758, -102.226526766947, 
    -102.221933499182, -102.217277175214, -102.21257059579, -102.207826701155, 
    -102.203058535141, -102.198279208964, -102.193501864824, 
    22.0859192364704, 22.085534238459, 22.0849172382595, 22.0840699348155, 
    22.0829946611303, 22.0816943777556, 22.0801726645306, 22.0784337105961, 
    22.0764823027115, 22.0743238119113, 22.0719641785371, 22.069409895692, 
    22.0666679911626, 22.0637460078641, 22.060651982863, 22.05739442504, 
    22.0539822914556, 22.0504249624898, 22.0467322158238, 22.0429141993409, 
    22.0389814030237, 22.034944629927, 22.0308149663093, 22.0266037510087, 
    22.0223225441472, 22.0179830952551, 22.0135973109014, 22.0091772219229, 
    22.0047349503427, 22.00028267607, 21.9958326034732, 21.9913969279179, 
    21.9869878023627, 21.9826173041033, 21.9782974017556, 21.974039922567, 
    21.9698565201456, 21.9657586426917, 21.9617575018195, 21.9578640420503, 
    21.9540889110588, 21.950442430752, 21.9469345692567, 21.9435749138888, 
    21.9403726451761, 21.9373365120024, 21.9344748079369, 21.9317953488116, 
    21.9293054516035, 21.9270119146773, 21.9249209994388, 21.9230384134463, 
    21.9213692950234, 21.9199181994126, 21.918689086505, 21.9176853101771, 
    21.9169096092626, 21.9163641001814, 21.9160502712448, 21.9159689786515, 
    21.9161204441842, 21.9165042546133, 21.9171193628082, 21.9179640905543, 
    21.9190361330684, 21.9203325652003, 21.9218498493057, 21.9235838447705, 
    21.9255298191611, 21.9276824609738, 21.9300358939496, 21.9325836929185, 
    21.9353189011326, 21.9382340490429, 21.9413211744729, 21.9445718441351, 
    21.947977176437, 21.951527865515, 21.9552142064358, 21.9590261214987, 
    21.9629531875695, 21.9669846643751, 21.9711095236842, 21.9753164792958, 
    21.9795940177575, 21.983930429731, 21.9883138419203, 21.9927322494785, 
    21.9971735488044, 22.001625570641, 22.0060761133856, 22.0105129765217, 
    22.0149239940793, 22.0192970680334, 22.0236202015486, 22.0278815319759, 
    22.0320693635126, 22.0361721994319, 22.0401787737937, 22.0440780825465, 
    22.0478594139329, 22.051512378113, 22.0550269359216, 22.0583934266765, 
    22.0616025949599, 22.0646456162951, 22.0675141216445, 22.0702002206599, 
    22.0726965236163, 22.0749961619661, 22.0770928074551, 22.0789806897431, 
    22.080654612479, 22.0821099677832, 22.0833427490944, 22.0843495623444, 
    22.0851276354261, 22.0856748259297, 22.0859896271218, 22.0860711721515, 
    22.0859192364704), .Dim = c(121L, 2L))), class = c("XY", 
    "POLYGON", "sfg")), structure(list(structure(c(-102.366703187571, 
    -102.361953401998, -102.357231241862, -102.352549690896, 
    -102.347921620104, -102.343359752058, -102.338876625611, 
    -102.334484561146, -102.330195626451, -102.326021603316, 
    -102.321973954954, -102.318063794337, -102.31430185353, -102.310698454116, 
    -102.307263478792, -102.304006344217, -102.300935975184, 
    -102.29806078019, -102.295388628469, -102.29292682855, -102.290682108398, 
    -102.288660597189, -102.286867808771, -102.285308626848, 
    -102.283987291923, -102.282907390044, -102.282071843364, 
    -102.281482902546, -102.281142141029, -102.281050451163, 
    -102.28120804222, -102.28161444028, -102.282268489991, -102.283168358188, 
    -102.284311539362, -102.285694862952, -102.287314502445, 
    -102.28916598624, -102.291244210262, -102.293543452269, -102.296057387819, 
    -102.298779107849, -102.301701137812, -102.304815458326, 
    -102.308113527261, -102.311586303219, -102.315224270331, 
    -102.319017464311, -102.322955499685, -102.327027598133, 
    -102.33122261786, -102.335529083924, -102.339935219426, -102.344428977502, 
    -102.348998074013, -102.353630020852, -102.358312159789, 
    -102.363031696755, -102.367775736476, -102.372531317382, 
    -102.37728544667, -102.38202513547, -102.386737433981, -102.391409466519, 
    -102.396028466363, -102.400581810326, -102.405057052946, 
    -102.409441960217, -102.413724542767, -102.417893088403, 
    -102.421936193923, -102.425842796126, -102.429602201935, 
    -102.43320411754, -102.4366386765, -102.439896466719, -102.442968556219, 
    -102.445846517653, -102.448522451477, -102.450989007725, 
    -102.453239406319, -102.45526745587, -102.457067570893, -102.458634787414, 
    -102.4599647769, -102.461053858482, -102.461899009434, -102.462497873866, 
    -102.462848769613, -102.462950693293, -102.462803323512, 
    -102.462407022203, -102.461762834097, -102.460872484312, 
    -102.459738374072, -102.458363574556, -102.456751818886, 
    -102.454907492285, -102.452835620412, -102.450541855915, 
    -102.448032463222, -102.445314301624, -102.442394806679, 
    -102.439281969997, -102.435984317449, -102.432510885871, 
    -102.428871198315, -102.425075237915, -102.421133420455, 
    -102.417056565692, -102.412855867531, -102.40854286313, -102.404129401022, 
    -102.399627608343, -102.395049857261, -102.390408730704, 
    -102.385716987477, -102.380987526877, -102.376233352903, 
    -102.371467538159, -102.366703187571, 21.8658255348988, 21.8654500907327, 
    21.8648428753359, 21.8640055607182, 21.8629404523928, 21.8616504829416, 
    21.8601392038325, 21.8584107755117, 21.8564699558012, 21.8543220866349, 
    21.8519730791717, 21.8494293973302, 21.8466980397919, 21.843786520526, 
    21.8407028478924, 21.8374555023831, 21.8340534130662, 21.8305059328018, 
    21.8268228122996, 21.8230141730931, 21.8190904795087, 21.8150625097083, 
    21.8109413258889, 21.8067382437227, 21.8024648011255, 21.7981327264397, 
    21.7937539061222, 21.7893403520278, 21.7849041683784, 21.7804575185102, 
    21.7760125914911, 21.7715815686993, 21.7671765904566, 21.7628097228052, 
    21.7584929245207, 21.7542380144494, 21.7500566392588, 21.7459602416876, 
    21.7419600293812, 21.7380669443942, 21.7342916334436, 21.7306444189888, 
    21.7271352712172, 21.7237737810075, 21.7205691339423, 21.7175300854374, 
    21.7146649370543, 21.7119815140545, 21.7094871442571, 21.7071886382512, 
    21.7050922710163, 21.7032037649975, 21.7015282746776, 21.7000703726872, 
    21.698834037488, 21.697822642659, 21.697038947816, 21.6964850911841, 
    21.6961625838445, 21.6960723056696, 21.6962145029558, 21.6965887877606, 
    21.6971941389471, 21.69802890493, 21.6990908081196, 21.7003769510501, 
    21.7018838241781, 21.70360731533, 21.7055427207767, 21.7076847579041, 
    21.7100275794504, 21.7125647892714, 21.7152894595944, 21.7181941497166, 
    21.7212709260998, 21.7245113838088, 21.7279066692395, 21.7314475040767, 
    21.7351242104195, 21.7389267370089, 21.7428446864898, 21.7468673436336, 
    21.7509837044501, 21.7551825061087, 21.7594522575916, 21.7637812709972, 
    21.7681576934083, 21.7725695392421, 21.7770047229928, 21.7814510922792, 
    21.7858964611076, 21.790328643259, 21.7947354857098, 21.799104901993, 
    21.8034249054097, 21.8076836419976, 21.811869423165, 21.8159707579007, 
    21.8199763844667, 21.8238753014886, 21.8276567983526, 21.8313104848251, 
    21.8348263198099, 21.8381946391619, 21.8414061824772, 21.844452118783, 
    21.8473240710538, 21.8500141394829, 21.8525149234429, 21.8548195420695, 
    21.8569216534119, 21.8588154720916, 21.8604957854188, 21.861957967921, 
    21.8631979942385, 21.8642124503512, 21.8649985431035, 21.8655541079981, 
    21.8658776152377, 21.8659681739947, 21.8658255348988), .Dim = c(121L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON", 
"sfc"), precision = 0, bbox = structure(c(xmin = -102.462950693293, 
ymin = 21.6960723056696, xmax = -102.107797978126, ymax = 22.0860711721515
), class = "bbox"), crs = structure(list(epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(station_number = NA_integer_, 
station_alt = NA_integer_, month = NA_integer_, Mean_min = NA_integer_, 
Mean_max = NA_integer_, months_observed = NA_integer_), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = 1:2, class = c("sf", 
"data.frame"))

I want to write a for loop that assigns the station_number to all the the polygons touched by buffer number i, then i + 1, etc... i.e. Assign station number 100 to all the polygons touched by staion 100's buffer.

To do this, I tried:

mex$station_number <- 0
for (i in 1:300) {
  mex[buffers[i, ], ]$station_number <- buffers[i, ]$station_number
}

However, running the loop (or even a named case) returns:

Error in Ops.sfc(left, right) : operation >= not supported

I expected this to do the following: The polygons touched by the ith buffer will be assigned the ith's station number.

I read that it was a problem with the $ in another forum. I tried escaping it using \ and \\ but it didn't work.

Arturo Sbr
  • 5,567
  • 4
  • 38
  • 76
  • It's hard to help without a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). There are example datasets that ship with `sf` that you could use here – camille Feb 19 '19 at 16:51
  • @camille thanks. I updated the question to display my actual data. I didn't do so initially because I thought it was a matter of syntax in my loop. The example should now be reproducible. You only have to set the CRS of the map from 5332 to 4326 using `st_set_crs(mex, 5332)` and `st_transform(mex, 4326)`. – Arturo Sbr Feb 19 '19 at 17:31
  • 1
    Try `st_join(mex, buffers)` when they have the same projection. – TimSalabim Feb 19 '19 at 18:48
  • @TimSalabim that did the trick. Thanks! – Arturo Sbr Feb 19 '19 at 23:07

1 Answers1

1

As pointed out by @TimSalabim, this can be done with st_join():

install.packages("lwgeom") # required when st_join(..., largest = TRUE)
mex <- st_join(x = mex, y = buffers, left = FALSE, largest = TRUE)

In fact, this solution is much more accurate than the loop I was writing because of the argument largest. The documentations states:

largest

logical; if TRUE, return x features augmented with the fields of y that have the largest overlap with each of the features of x

So basically, this leaves me with the buffer that has the largest intersection over a given polygon.

Community
  • 1
  • 1
Arturo Sbr
  • 5,567
  • 4
  • 38
  • 76