Spatial Query (Selection) with Apply

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Spatial Query (Selection) with Apply

aestfd
Hi,

I want to do multiple selections of a point shapefile based on polygons on
other layers, I can do this in a for loop, but I desire to do this in a
function of the apply family.

I named the point shapefile "nodes" and the polygons shapefile "zones".

This is what I did:

tmp <- list()
for (i in 1:nrow(zones@data)) {
    tmp[[i]] <- nodes[subset(zones, ESTUDIO == i),]
    tmp
  }

But I have no clue how to change it to the apply family, can you provide
an example of this?

Thanks in advance.

Regards,
Ariel Fuentes

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Spatial Query (Selection) with Apply

robinlovelace
Hi Ariel,

It helps when asking for help with code to produce a reproducible example.
To understand your input data I've created nodes and zones based on the
spData data nz_height and nz_elev. Based on the assumption you want to find
all the nodes in each zone I think the direct answer to your question is
something like the following:

tmp2 = lapply(s, function(x) {
  nodes[zones[x, ], ]
})

The longer answer is that aggregate(), st_join() + aggregate()/summarize()
may provide quicker solutions, depending on what you want to do with the
points after grouping them by which zone they fall in.
Note: I've used sf objects based on the explanation here https://geocompr.
robinlovelace.net/spatial-operations.html#spatial-vec which may not work
with sp data but the aggregate code should work roughly the same:

# install.packages("sf")
# install.packages("spData")
library(spData)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
zones = nz
nodes = nz_height
tmp = list()
s = 1:nrow(nz)

for(i in s) {
  tmp[[i]] = nodes[nz[i, ], ]
}

# understand what's going on with plots (not shown)
# plot(st_geometry(nz))
# plot(st_geometry(tmp[[i]]), add = TRUE, col = "red")
# plot(nz[i, ], col = "green", add = TRUE)

tmp2 = lapply(s, function(x) {
  nz_height[nz[x, ], ]
})

identical(tmp, tmp2)
#> [1] TRUE

Created on 2018-08-27 by the [reprex package](http://reprex.tidyverse.org)
(v0.2.0).

I've also pasted the reprex into the geocompr github tracker so the plots
can be seen and the code formatted: https://github.com/
Robinlovelace/geocompr/issues/294

Hope this helps,

Robin


On Mon, Aug 27, 2018 at 9:53 PM, Ariel Fuentesdi <[hidden email]>
wrote:

> Hi,
>
> I want to do multiple selections of a point shapefile based on polygons on
> other layers, I can do this in a for loop, but I desire to do this in a
> function of the apply family.
>
> I named the point shapefile "nodes" and the polygons shapefile "zones".
>
> This is what I did:
>
> tmp <- list()
> for (i in 1:nrow(zones@data)) {
>     tmp[[i]] <- nodes[subset(zones, ESTUDIO == i),]
>     tmp
>   }
>
> But I have no clue how to change it to the apply family, can you provide
> an example of this?
>
> Thanks in advance.
>
> Regards,
> Ariel Fuentes
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Spatial Query (Selection) with Apply

aestfd
Thanks, Robin

It worked perfectly, although I transformed the data from sf to
SpatialDataFrame objects with:

zones <- as(zones, "Spatial")
nodes <- as(nodes, "Spatial")

The next step I will follow is to create Distance Matrix like this:

coord <- tmp2[[1]]@coords
dist1 <- apply(X = coord, MARGIN = 1, FUN = function(x) spDistsN1(coord, x,
longlat = T))

obviously, the task is a distance matrix for every selection. ¿What would
you recommend, the aggregate function as you talked before or something
like a nested apply?

Regards,
Ariel



2018-08-27 19:32 GMT-03:00 Robin Lovelace <[hidden email]>:

> Hi Ariel,
>
> It helps when asking for help with code to produce a reproducible example.
> To understand your input data I've created nodes and zones based on the
> spData data nz_height and nz_elev. Based on the assumption you want to find
> all the nodes in each zone I think the direct answer to your question is
> something like the following:
>
> tmp2 = lapply(s, function(x) {
>   nodes[zones[x, ], ]
> })
>
> The longer answer is that aggregate(), st_join() + aggregate()/summarize()
> may provide quicker solutions, depending on what you want to do with the
> points after grouping them by which zone they fall in.
> Note: I've used sf objects based on the explanation here
> https://geocompr.robinlovelace.net/spatial-operations.html#spatial-vec
> which may not work with sp data but the aggregate code should work roughly
> the same:
>
> # install.packages("sf")
> # install.packages("spData")
> library(spData)
> library(sf)
> #> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
> zones = nz
> nodes = nz_height
> tmp = list()
> s = 1:nrow(nz)
>
> for(i in s) {
>   tmp[[i]] = nodes[nz[i, ], ]
> }
>
> # understand what's going on with plots (not shown)
> # plot(st_geometry(nz))
> # plot(st_geometry(tmp[[i]]), add = TRUE, col = "red")
> # plot(nz[i, ], col = "green", add = TRUE)
>
> tmp2 = lapply(s, function(x) {
>   nz_height[nz[x, ], ]
> })
>
> identical(tmp, tmp2)
> #> [1] TRUE
>
> Created on 2018-08-27 by the [reprex package](http://reprex.tidyverse.org)
> (v0.2.0).
>
> I've also pasted the reprex into the geocompr github tracker so the plots
> can be seen and the code formatted: https://github.com/Robinlovela
> ce/geocompr/issues/294
>
> Hope this helps,
>
> Robin
>
>
> On Mon, Aug 27, 2018 at 9:53 PM, Ariel Fuentesdi <[hidden email]
> > wrote:
>
>> Hi,
>>
>> I want to do multiple selections of a point shapefile based on polygons on
>> other layers, I can do this in a for loop, but I desire to do this in a
>> function of the apply family.
>>
>> I named the point shapefile "nodes" and the polygons shapefile "zones".
>>
>> This is what I did:
>>
>> tmp <- list()
>> for (i in 1:nrow(zones@data)) {
>>     tmp[[i]] <- nodes[subset(zones, ESTUDIO == i),]
>>     tmp
>>   }
>>
>> But I have no clue how to change it to the apply family, can you provide
>> an example of this?
>>
>> Thanks in advance.
>>
>> Regards,
>> Ariel Fuentes
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Spatial Query (Selection) with Apply

robinlovelace
Hi Ariel,

Glad it helped.

Yes you could place your apply command into a lapply to get a list of
matrices.

Something like this should work:

dist_list <- lapply(tmp, function(coord) {

   apply(X = coord, MARGIN = 1, FUN = function(x) spDistsN1(coord, x,
longlat = T)

}

Robin



On Tue, Aug 28, 2018 at 9:13 PM, Ariel Fuentesdi <[hidden email]>
wrote:

> Thanks, Robin
>
> It worked perfectly, although I transformed the data from sf to
> SpatialDataFrame objects with:
>
> zones <- as(zones, "Spatial")
> nodes <- as(nodes, "Spatial")
>
> The next step I will follow is to create Distance Matrix like this:
>
> coord <- tmp2[[1]]@coords
> dist1 <- apply(X = coord, MARGIN = 1, FUN = function(x) spDistsN1(coord,
> x, longlat = T))
>
> obviously, the task is a distance matrix for every selection. ¿What would
> you recommend, the aggregate function as you talked before or something
> like a nested apply?
>
> Regards,
> Ariel
>
>
>
> 2018-08-27 19:32 GMT-03:00 Robin Lovelace <[hidden email]>:
>
>> Hi Ariel,
>>
>> It helps when asking for help with code to produce a reproducible
>> example. To understand your input data I've created nodes and zones based
>> on the spData data nz_height and nz_elev. Based on the assumption you want
>> to find all the nodes in each zone I think the direct answer to your
>> question is something like the following:
>>
>> tmp2 = lapply(s, function(x) {
>>   nodes[zones[x, ], ]
>> })
>>
>> The longer answer is that aggregate(), st_join() +
>> aggregate()/summarize() may provide quicker solutions, depending on what
>> you want to do with the points after grouping them by which zone they fall
>> in.
>> Note: I've used sf objects based on the explanation here
>> https://geocompr.robinlovelace.net/spatial-operations.html#spatial-vec
>> which may not work with sp data but the aggregate code should work roughly
>> the same:
>>
>> # install.packages("sf")
>> # install.packages("spData")
>> library(spData)
>> library(sf)
>> #> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
>> zones = nz
>> nodes = nz_height
>> tmp = list()
>> s = 1:nrow(nz)
>>
>> for(i in s) {
>>   tmp[[i]] = nodes[nz[i, ], ]
>> }
>>
>> # understand what's going on with plots (not shown)
>> # plot(st_geometry(nz))
>> # plot(st_geometry(tmp[[i]]), add = TRUE, col = "red")
>> # plot(nz[i, ], col = "green", add = TRUE)
>>
>> tmp2 = lapply(s, function(x) {
>>   nz_height[nz[x, ], ]
>> })
>>
>> identical(tmp, tmp2)
>> #> [1] TRUE
>>
>> Created on 2018-08-27 by the [reprex package](http://reprex.tidyverse.org)
>> (v0.2.0).
>>
>> I've also pasted the reprex into the geocompr github tracker so the plots
>> can be seen and the code formatted: https://github.com/Robinlovela
>> ce/geocompr/issues/294
>>
>> Hope this helps,
>>
>> Robin
>>
>>
>> On Mon, Aug 27, 2018 at 9:53 PM, Ariel Fuentesdi <
>> [hidden email]> wrote:
>>
>>> Hi,
>>>
>>> I want to do multiple selections of a point shapefile based on polygons
>>> on
>>> other layers, I can do this in a for loop, but I desire to do this in a
>>> function of the apply family.
>>>
>>> I named the point shapefile "nodes" and the polygons shapefile "zones".
>>>
>>> This is what I did:
>>>
>>> tmp <- list()
>>> for (i in 1:nrow(zones@data)) {
>>>     tmp[[i]] <- nodes[subset(zones, ESTUDIO == i),]
>>>     tmp
>>>   }
>>>
>>> But I have no clue how to change it to the apply family, can you provide
>>> an example of this?
>>>
>>> Thanks in advance.
>>>
>>> Regards,
>>> Ariel Fuentes
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo