This is an

**for**__archive__**R-sig-geo**, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first**(through the list homepage, not via Nabble)**and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.*Updated:*1 hour 57 min ago

### How to draw a legend when plotting sf (simple feature) objects

Dear list members,

How can I draw a legend when plotting an sf (simple feature) object.

nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)

plot(nc["SID79"])

Best,

--

*Manuel Spínola, Ph.D.*

Instituto Internacional en Conservación y Manejo de Vida Silvestre

Universidad Nacional

Apartado 1350-3000

Heredia

COSTA RICA

[hidden email] <[hidden email]>

[hidden email]

Teléfono: (506) 8706 - 4662

Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>

Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

How can I draw a legend when plotting an sf (simple feature) object.

nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)

plot(nc["SID79"])

Best,

--

*Manuel Spínola, Ph.D.*

Instituto Internacional en Conservación y Manejo de Vida Silvestre

Universidad Nacional

Apartado 1350-3000

Heredia

COSTA RICA

[hidden email] <[hidden email]>

[hidden email]

Teléfono: (506) 8706 - 4662

Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>

Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: rasterTopolygons: singleparts often needed

Not at my computer, but does sp::disaggregate(v) return your singlepart

polygons?

On Sun, 26 Feb 2017 at 22:40, Agustin Lobo <[hidden email]> wrote:

> Could raster::rasterTopolygons() be modified in future versions

>

> to optionally output singlepart polygons? e.g, in the following

>

> r <-raster(nrow=16, ncol=16)

>

> d <- rep(rep(1:4,rep(4,4)),8)

>

> r[] <- c(d,rev(d))

>

> plot(r)

>

> v <- rasterToPolygons(r,dissolve=TRUE)

>

> plot(v[1,])

>

> dim(v@data)

>

>

>

> the new option

>

> v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)

>

> would result on

>

> plot(v[1,])

>

> displaying one single square and

>

> dim(v@data)

>

> would be 8 rows and 1 column

>

>

>

> (equivalent to gdal_polygonize)

>

>

>

> Agus

>

>

>

> _______________________________________________

>

> 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

polygons?

On Sun, 26 Feb 2017 at 22:40, Agustin Lobo <[hidden email]> wrote:

> Could raster::rasterTopolygons() be modified in future versions

>

> to optionally output singlepart polygons? e.g, in the following

>

> r <-raster(nrow=16, ncol=16)

>

> d <- rep(rep(1:4,rep(4,4)),8)

>

> r[] <- c(d,rev(d))

>

> plot(r)

>

> v <- rasterToPolygons(r,dissolve=TRUE)

>

> plot(v[1,])

>

> dim(v@data)

>

>

>

> the new option

>

> v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)

>

> would result on

>

> plot(v[1,])

>

> displaying one single square and

>

> dim(v@data)

>

> would be 8 rows and 1 column

>

>

>

> (equivalent to gdal_polygonize)

>

>

>

> Agus

>

>

>

> _______________________________________________

>

> 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

### Re: Assign a unique ID number to each patch in a raster

Yes, my gist mentioned by Agustin should do the trick. See the example in

the comment at https://gist.github.com/johnbaums/6d155cfca28e02b05ad5.

On Sun, 26 Feb 2017 at 21:42, Agustin Lobo <[hidden email]> wrote:

> I think Nelly points out that package raster lacks a more general

>

> clumping tool (i.e., the equivalent to r.clump in grass) that I've

>

> often missed too.

>

> I'm not sure circumventing the way Ben suggest would work for large

>

> raster layers, for which package raster was designed.

>

> I have a note to myself regarding this

>

> https://gist.github.com/johnbaums/6d155cfca28e02b05ad5

>

> and

>

> www.openforis.org/tools/geospatial-toolkit.html

>

> but have not tried it.

>

>

>

> Agus

>

>

>

> On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:

>

> > Hi,

>

> >

>

> > I think you are asking for connected component labeling. raster::clump()

> expects your features to be like islands surrounded by NA or 0 background.

> There might be more Raster* centric ways to do it, but you could use plain

> old vanilla image processing using the very good imager package (

> https://cran.r-project.org/web/packages/imager/ <

> https://cran.r-project.org/web/packages/imager/> ). The steps below show

> in color your actual values with the connected component IDs labeled on

> each cell. I used the 8-connected labeling but 4-connected is also

> available.

>

> >

>

> > library(raster)

>

> > library(imager)

>

> >

>

> > set.seed(10)

>

> > nc <- 8

>

> > nr <- 8

>

> > MIN <- 1

>

> > MAX <- 4

>

> > LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")

>

> >

>

> > r <- raster(ncols=nc, nrows=nr)

>

> > r[] <- round(runif(ncell(r),MIN, MAX),digits=0)

>

> > cc <- r

>

> >

>

> > img <- imager::as.cimg(as.matrix(r))

>

> > labeled <- imager::label(img, high_connectivity = TRUE)

>

> > cc[] <- as.matrix(labeled)

>

> >

>

> > plot(r, col = LUT)

>

> > txt <- cc[]

>

> > xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))

>

> > text(xy[,1], xy[,2], txt )

>

> >

>

> >

>

> > Cheers,

>

> > Ben

>

> >> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:

>

> >>

>

> >> Hello,

>

> >> I would like to assign a unique ID number to each patch (a patch is

> composed of a set of adjacent cells) as shown in this figure:

>

> >> <2328.png>

>

> >>

>

> >> I tested the clump function (package raster) by applying an eight

> adjacent cells rule but this function mixes all cell values into single

> patch.

>

> >>

>

> >> Here is a code example to create a raster. The raster contains values

> ranged from 1 to 9.

>

> >> r <- raster(ncols=12, nrows=12)

>

> >> r[] <- round(runif(ncell(r),1,9),digits=0)

>

> >> plot(r)

>

> >>

>

> >> Is there a way to assign a unique ID number to grouping cells that form

> a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?

>

> >>

>

> >> Thanks a lot for your time

>

> >> Nell

>

> >>

>

> >>

>

> >> _______________________________________________

>

> >> R-sig-Geo mailing list

>

> >> [hidden email] <mailto:[hidden email]>

>

> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo>

>

> > Ben Tupper

>

> > Bigelow Laboratory for Ocean Sciences

>

> > 60 Bigelow Drive, P.O. Box 380

>

> > East Boothbay, Maine 04544

>

> > http://www.bigelow.org

>

> >

>

> >

>

> >

>

> >

>

> > [[alternative HTML version deleted]]

>

> >

>

> > _______________________________________________

>

> > R-sig-Geo mailing list

>

> > [hidden email]

>

> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

>

>

> _______________________________________________

>

> 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

the comment at https://gist.github.com/johnbaums/6d155cfca28e02b05ad5.

On Sun, 26 Feb 2017 at 21:42, Agustin Lobo <[hidden email]> wrote:

> I think Nelly points out that package raster lacks a more general

>

> clumping tool (i.e., the equivalent to r.clump in grass) that I've

>

> often missed too.

>

> I'm not sure circumventing the way Ben suggest would work for large

>

> raster layers, for which package raster was designed.

>

> I have a note to myself regarding this

>

> https://gist.github.com/johnbaums/6d155cfca28e02b05ad5

>

> and

>

> www.openforis.org/tools/geospatial-toolkit.html

>

> but have not tried it.

>

>

>

> Agus

>

>

>

> On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:

>

> > Hi,

>

> >

>

> > I think you are asking for connected component labeling. raster::clump()

> expects your features to be like islands surrounded by NA or 0 background.

> There might be more Raster* centric ways to do it, but you could use plain

> old vanilla image processing using the very good imager package (

> https://cran.r-project.org/web/packages/imager/ <

> https://cran.r-project.org/web/packages/imager/> ). The steps below show

> in color your actual values with the connected component IDs labeled on

> each cell. I used the 8-connected labeling but 4-connected is also

> available.

>

> >

>

> > library(raster)

>

> > library(imager)

>

> >

>

> > set.seed(10)

>

> > nc <- 8

>

> > nr <- 8

>

> > MIN <- 1

>

> > MAX <- 4

>

> > LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")

>

> >

>

> > r <- raster(ncols=nc, nrows=nr)

>

> > r[] <- round(runif(ncell(r),MIN, MAX),digits=0)

>

> > cc <- r

>

> >

>

> > img <- imager::as.cimg(as.matrix(r))

>

> > labeled <- imager::label(img, high_connectivity = TRUE)

>

> > cc[] <- as.matrix(labeled)

>

> >

>

> > plot(r, col = LUT)

>

> > txt <- cc[]

>

> > xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))

>

> > text(xy[,1], xy[,2], txt )

>

> >

>

> >

>

> > Cheers,

>

> > Ben

>

> >> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:

>

> >>

>

> >> Hello,

>

> >> I would like to assign a unique ID number to each patch (a patch is

> composed of a set of adjacent cells) as shown in this figure:

>

> >> <2328.png>

>

> >>

>

> >> I tested the clump function (package raster) by applying an eight

> adjacent cells rule but this function mixes all cell values into single

> patch.

>

> >>

>

> >> Here is a code example to create a raster. The raster contains values

> ranged from 1 to 9.

>

> >> r <- raster(ncols=12, nrows=12)

>

> >> r[] <- round(runif(ncell(r),1,9),digits=0)

>

> >> plot(r)

>

> >>

>

> >> Is there a way to assign a unique ID number to grouping cells that form

> a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?

>

> >>

>

> >> Thanks a lot for your time

>

> >> Nell

>

> >>

>

> >>

>

> >> _______________________________________________

>

> >> R-sig-Geo mailing list

>

> >> [hidden email] <mailto:[hidden email]>

>

> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo>

>

> > Ben Tupper

>

> > Bigelow Laboratory for Ocean Sciences

>

> > 60 Bigelow Drive, P.O. Box 380

>

> > East Boothbay, Maine 04544

>

> > http://www.bigelow.org

>

> >

>

> >

>

> >

>

> >

>

> > [[alternative HTML version deleted]]

>

> >

>

> > _______________________________________________

>

> > R-sig-Geo mailing list

>

> > [hidden email]

>

> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

>

>

> _______________________________________________

>

> 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

### rasterTopolygons: singleparts often needed

Could raster::rasterTopolygons() be modified in future versions

to optionally output singlepart polygons? e.g, in the following

r <-raster(nrow=16, ncol=16)

d <- rep(rep(1:4,rep(4,4)),8)

r[] <- c(d,rev(d))

plot(r)

v <- rasterToPolygons(r,dissolve=TRUE)

plot(v[1,])

dim(v@data)

the new option

v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)

would result on

plot(v[1,])

displaying one single square and

dim(v@data)

would be 8 rows and 1 column

(equivalent to gdal_polygonize)

Agus

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

to optionally output singlepart polygons? e.g, in the following

r <-raster(nrow=16, ncol=16)

d <- rep(rep(1:4,rep(4,4)),8)

r[] <- c(d,rev(d))

plot(r)

v <- rasterToPolygons(r,dissolve=TRUE)

plot(v[1,])

dim(v@data)

the new option

v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)

would result on

plot(v[1,])

displaying one single square and

dim(v@data)

would be 8 rows and 1 column

(equivalent to gdal_polygonize)

Agus

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: Assign a unique ID number to each patch in a raster

I think Nelly points out that package raster lacks a more general

clumping tool (i.e., the equivalent to r.clump in grass) that I've

often missed too.

I'm not sure circumventing the way Ben suggest would work for large

raster layers, for which package raster was designed.

I have a note to myself regarding this

https://gist.github.com/johnbaums/6d155cfca28e02b05ad5

and

www.openforis.org/tools/geospatial-toolkit.html

but have not tried it.

Agus

On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:

> Hi,

>

> I think you are asking for connected component labeling. raster::clump() expects your features to be like islands surrounded by NA or 0 background. There might be more Raster* centric ways to do it, but you could use plain old vanilla image processing using the very good imager package ( https://cran.r-project.org/web/packages/imager/ <https://cran.r-project.org/web/packages/imager/> ). The steps below show in color your actual values with the connected component IDs labeled on each cell. I used the 8-connected labeling but 4-connected is also available.

>

> library(raster)

> library(imager)

>

> set.seed(10)

> nc <- 8

> nr <- 8

> MIN <- 1

> MAX <- 4

> LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")

>

> r <- raster(ncols=nc, nrows=nr)

> r[] <- round(runif(ncell(r),MIN, MAX),digits=0)

> cc <- r

>

> img <- imager::as.cimg(as.matrix(r))

> labeled <- imager::label(img, high_connectivity = TRUE)

> cc[] <- as.matrix(labeled)

>

> plot(r, col = LUT)

> txt <- cc[]

> xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))

> text(xy[,1], xy[,2], txt )

>

>

> Cheers,

> Ben

>> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:

>>

>> Hello,

>> I would like to assign a unique ID number to each patch (a patch is composed of a set of adjacent cells) as shown in this figure:

>> <2328.png>

>>

>> I tested the clump function (package raster) by applying an eight adjacent cells rule but this function mixes all cell values into single patch.

>>

>> Here is a code example to create a raster. The raster contains values ranged from 1 to 9.

>> r <- raster(ncols=12, nrows=12)

>> r[] <- round(runif(ncell(r),1,9),digits=0)

>> plot(r)

>>

>> Is there a way to assign a unique ID number to grouping cells that form a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?

>>

>> Thanks a lot for your time

>> Nell

>>

>>

>> _______________________________________________

>> R-sig-Geo mailing list

>> [hidden email] <mailto:[hidden email]>

>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>

> Ben Tupper

> Bigelow Laboratory for Ocean Sciences

> 60 Bigelow Drive, P.O. Box 380

> East Boothbay, Maine 04544

> http://www.bigelow.org

>

>

>

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

clumping tool (i.e., the equivalent to r.clump in grass) that I've

often missed too.

I'm not sure circumventing the way Ben suggest would work for large

raster layers, for which package raster was designed.

I have a note to myself regarding this

https://gist.github.com/johnbaums/6d155cfca28e02b05ad5

and

www.openforis.org/tools/geospatial-toolkit.html

but have not tried it.

Agus

On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:

> Hi,

>

> I think you are asking for connected component labeling. raster::clump() expects your features to be like islands surrounded by NA or 0 background. There might be more Raster* centric ways to do it, but you could use plain old vanilla image processing using the very good imager package ( https://cran.r-project.org/web/packages/imager/ <https://cran.r-project.org/web/packages/imager/> ). The steps below show in color your actual values with the connected component IDs labeled on each cell. I used the 8-connected labeling but 4-connected is also available.

>

> library(raster)

> library(imager)

>

> set.seed(10)

> nc <- 8

> nr <- 8

> MIN <- 1

> MAX <- 4

> LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")

>

> r <- raster(ncols=nc, nrows=nr)

> r[] <- round(runif(ncell(r),MIN, MAX),digits=0)

> cc <- r

>

> img <- imager::as.cimg(as.matrix(r))

> labeled <- imager::label(img, high_connectivity = TRUE)

> cc[] <- as.matrix(labeled)

>

> plot(r, col = LUT)

> txt <- cc[]

> xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))

> text(xy[,1], xy[,2], txt )

>

>

> Cheers,

> Ben

>> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:

>>

>> Hello,

>> I would like to assign a unique ID number to each patch (a patch is composed of a set of adjacent cells) as shown in this figure:

>> <2328.png>

>>

>> I tested the clump function (package raster) by applying an eight adjacent cells rule but this function mixes all cell values into single patch.

>>

>> Here is a code example to create a raster. The raster contains values ranged from 1 to 9.

>> r <- raster(ncols=12, nrows=12)

>> r[] <- round(runif(ncell(r),1,9),digits=0)

>> plot(r)

>>

>> Is there a way to assign a unique ID number to grouping cells that form a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?

>>

>> Thanks a lot for your time

>> Nell

>>

>>

>> _______________________________________________

>> R-sig-Geo mailing list

>> [hidden email] <mailto:[hidden email]>

>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>

> Ben Tupper

> Bigelow Laboratory for Ocean Sciences

> 60 Bigelow Drive, P.O. Box 380

> East Boothbay, Maine 04544

> http://www.bigelow.org

>

>

>

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### raster brick from netcdf - specify range

Want to read part of a netcdf file into a raster brick. The file (in

ncdf4 order) contains a variable

pre[lon, lat, time]

and those dimensions are: 720, 360, 1356

The following reads in all time layers:

rb <- brick( filename, varname=varname, nl=2 )

> rb

class : RasterBrick

dimensions : 360, 720, 259200, 1356 (nrow, ncol, ncell, nlayers)

resolution : 0.5, 0.5 (x, y)

extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

z-value : 12, 1367 (min, max)

varname : pre

Likewise, I've had no success limiting the read using xmn, xmx, ymn,

ymx, or template. The help page for brick suggests that all these should

work, or at least, it does not seem say that they don't work for netcdf

files.

My workaround requires directly reading the netcdf file, then making a

brick from that array, but when reading finer resolution data, the

memory use and speed could be an issue due to creating 2 objects instead

of 1.

raster version = 2.5-8

I'd appreciate any suggestions about this.

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

ncdf4 order) contains a variable

pre[lon, lat, time]

and those dimensions are: 720, 360, 1356

The following reads in all time layers:

rb <- brick( filename, varname=varname, nl=2 )

> rb

class : RasterBrick

dimensions : 360, 720, 259200, 1356 (nrow, ncol, ncell, nlayers)

resolution : 0.5, 0.5 (x, y)

extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

z-value : 12, 1367 (min, max)

varname : pre

Likewise, I've had no success limiting the read using xmn, xmx, ymn,

ymx, or template. The help page for brick suggests that all these should

work, or at least, it does not seem say that they don't work for netcdf

files.

My workaround requires directly reading the netcdf file, then making a

brick from that array, but when reading finer resolution data, the

memory use and speed could be an issue due to creating 2 objects instead

of 1.

raster version = 2.5-8

I'd appreciate any suggestions about this.

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: split SLDF with attributes

Marta,

Since you have intersecting lines in your SLDF, another approach might

be to use a loop (!) through your SpatialLines objects and carry the

attributes, something like:

library(raster)

sldf_in <- shapefile("mergedDF_1_DF_short")

# Init a list to keep all segments

sldf_list <- vector("list", nrow(sldf_in))

names(sldf_list) <- row.names(sldf_in)

for (i in names(sldf_list)) {

sl <- SegmentSpatialLines(as(sldf_in[i,], "SpatialLines"),

length=10000, merge.last=TRUE)

df <- sldf_in[i,]@data

df <- df[rep(seq_len(nrow(df)), length(sl)), ]

df$ID <- i

sldf_list[[i]] <- SpatialLinesDataFrame(sl, df, match.ID=FALSE)

# Make arbitrary unique IDs (for combining in next step)

sldf_list[[i]] <- spChFIDs(sldf_list[[i]], paste(i,

row.names(sldf_list[[i]]), sep="-"))

}

# Combine into a unique SLDF

sldf_out <- do.call(rbind, sldf_list)

--Mel.

On 2/24/2017 10:31 AM, marta azores wrote:

> Thanks for your answer Mel,

> Option1) Createsegments: this script cut the SLDF...but as autoridades

> said um the website...The results are spatiallines, not SLDF. Then,

> they don't jeep the attributes.

> Option2) rgeos:over I know this tool, but O have several lines

> crossing among them. I am afraid, that can be a bit messy.

> I doesn't know the gTouches...I check it after your suggestion..

> But I think the rgeos:Over os the best option for me.

> Thanks for your suggestion

> Marta

>

>

>

> No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie

> <[hidden email] <mailto:[hidden email]>> escreveu:

>

> Marta,

> Option 1) using CreateSegments() seems to do what you want. Why

> don't you simply use rgeos::over() or rgeos::gTouches() between

> your old and new SLDFs to match the index positions of your old

> attributes to your new segments?

>

> --Mel.

>

>

> On 2/23/2017 2:27 PM, marta azores wrote:

>

> Hi all,

>

> I am working with a large datasets of boat positions, and

> based on this

> points I created a SpatialLinesDataFrame. Now, I would like to

> split/divide SpatialLinesDataframe

> into n segments of 10km and keep the attributes. So far, I

> have tried

> multiple things:

>

> 1. The function from this website:

> http://rstudio-pubs-static.s3.

> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html

> <http://amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html>.

> The problem with

> it is that the function cuts the lines in segments of 10km but

> it is not

> possible to get the attributes.

>

> 2. I have also tried the function from this website:

> http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

> <http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline->

> shapefiles-into-equally-length-smaller-segments. The problem

> is that do not

> cut the segments in 10km, and it seems to cut the segments

> with random

> lengths.

>

> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

> <http://nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#>

> a7583629 propose to divide lines by segments but do not

> specify how it is

> possible to do in based on the length of the segment.

>

> Is there any function to do this in R? Or does somebody know

> the solution

> for that?

>

> My data looks like this:

>

> #Script:

>

> path_data <- "C:/R/"

>

> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>

>

> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

> <https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing>

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

> <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

Since you have intersecting lines in your SLDF, another approach might

be to use a loop (!) through your SpatialLines objects and carry the

attributes, something like:

library(raster)

sldf_in <- shapefile("mergedDF_1_DF_short")

# Init a list to keep all segments

sldf_list <- vector("list", nrow(sldf_in))

names(sldf_list) <- row.names(sldf_in)

for (i in names(sldf_list)) {

sl <- SegmentSpatialLines(as(sldf_in[i,], "SpatialLines"),

length=10000, merge.last=TRUE)

df <- sldf_in[i,]@data

df <- df[rep(seq_len(nrow(df)), length(sl)), ]

df$ID <- i

sldf_list[[i]] <- SpatialLinesDataFrame(sl, df, match.ID=FALSE)

# Make arbitrary unique IDs (for combining in next step)

sldf_list[[i]] <- spChFIDs(sldf_list[[i]], paste(i,

row.names(sldf_list[[i]]), sep="-"))

}

# Combine into a unique SLDF

sldf_out <- do.call(rbind, sldf_list)

--Mel.

On 2/24/2017 10:31 AM, marta azores wrote:

> Thanks for your answer Mel,

> Option1) Createsegments: this script cut the SLDF...but as autoridades

> said um the website...The results are spatiallines, not SLDF. Then,

> they don't jeep the attributes.

> Option2) rgeos:over I know this tool, but O have several lines

> crossing among them. I am afraid, that can be a bit messy.

> I doesn't know the gTouches...I check it after your suggestion..

> But I think the rgeos:Over os the best option for me.

> Thanks for your suggestion

> Marta

>

>

>

> No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie

> <[hidden email] <mailto:[hidden email]>> escreveu:

>

> Marta,

> Option 1) using CreateSegments() seems to do what you want. Why

> don't you simply use rgeos::over() or rgeos::gTouches() between

> your old and new SLDFs to match the index positions of your old

> attributes to your new segments?

>

> --Mel.

>

>

> On 2/23/2017 2:27 PM, marta azores wrote:

>

> Hi all,

>

> I am working with a large datasets of boat positions, and

> based on this

> points I created a SpatialLinesDataFrame. Now, I would like to

> split/divide SpatialLinesDataframe

> into n segments of 10km and keep the attributes. So far, I

> have tried

> multiple things:

>

> 1. The function from this website:

> http://rstudio-pubs-static.s3.

> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html

> <http://amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html>.

> The problem with

> it is that the function cuts the lines in segments of 10km but

> it is not

> possible to get the attributes.

>

> 2. I have also tried the function from this website:

> http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

> <http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline->

> shapefiles-into-equally-length-smaller-segments. The problem

> is that do not

> cut the segments in 10km, and it seems to cut the segments

> with random

> lengths.

>

> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

> <http://nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#>

> a7583629 propose to divide lines by segments but do not

> specify how it is

> possible to do in based on the length of the segment.

>

> Is there any function to do this in R? Or does somebody know

> the solution

> for that?

>

> My data looks like this:

>

> #Script:

>

> path_data <- "C:/R/"

>

> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>

>

> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

> <https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing>

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

> <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

### Re: split SLDF with attributes

Thanks for your answer Mel,

Option1) Createsegments: this script cut the SLDF...but as autoridades said

um the website...The results are spatiallines, not SLDF. Then, they don't

jeep the attributes.

Option2) rgeos:over I know this tool, but O have several lines crossing

among them. I am afraid, that can be a bit messy.

I doesn't know the gTouches...I check it after your suggestion..

But I think the rgeos:Over os the best option for me.

Thanks for your suggestion

Marta

No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie <[hidden email]>

escreveu:

> Marta,

> Option 1) using CreateSegments() seems to do what you want. Why don't you

> simply use rgeos::over() or rgeos::gTouches() between your old and new

> SLDFs to match the index positions of your old attributes to your new

> segments?

>

> --Mel.

>

>

> On 2/23/2017 2:27 PM, marta azores wrote:

>

>> Hi all,

>>

>> I am working with a large datasets of boat positions, and based on this

>> points I created a SpatialLinesDataFrame. Now, I would like to

>> split/divide SpatialLinesDataframe

>> into n segments of 10km and keep the attributes. So far, I have tried

>> multiple things:

>>

>> 1. The function from this website: http://rstudio-pubs-static.s3.

>> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem

>> with

>> it is that the function cuts the lines in segments of 10km but it is not

>> possible to get the attributes.

>>

>> 2. I have also tried the function from this website:

>> http://stackoverflow.com/questions/38700246/how-do-i-split-

>> divide-polyline-

>> shapefiles-into-equally-length-smaller-segments. The problem is that do

>> not

>> cut the segments in 10km, and it seems to cut the segments with random

>> lengths.

>>

>> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

>> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

>> a7583629 propose to divide lines by segments but do not specify how it is

>> possible to do in based on the length of the segment.

>>

>> Is there any function to do this in R? Or does somebody know the solution

>> for that?

>>

>> My data looks like this:

>>

>> #Script:

>>

>> path_data <- "C:/R/"

>>

>> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>>

>>

>> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklp

>> d3ZnRzA?usp=sharing

>>

>> [[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

Option1) Createsegments: this script cut the SLDF...but as autoridades said

um the website...The results are spatiallines, not SLDF. Then, they don't

jeep the attributes.

Option2) rgeos:over I know this tool, but O have several lines crossing

among them. I am afraid, that can be a bit messy.

I doesn't know the gTouches...I check it after your suggestion..

But I think the rgeos:Over os the best option for me.

Thanks for your suggestion

Marta

No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie <[hidden email]>

escreveu:

> Marta,

> Option 1) using CreateSegments() seems to do what you want. Why don't you

> simply use rgeos::over() or rgeos::gTouches() between your old and new

> SLDFs to match the index positions of your old attributes to your new

> segments?

>

> --Mel.

>

>

> On 2/23/2017 2:27 PM, marta azores wrote:

>

>> Hi all,

>>

>> I am working with a large datasets of boat positions, and based on this

>> points I created a SpatialLinesDataFrame. Now, I would like to

>> split/divide SpatialLinesDataframe

>> into n segments of 10km and keep the attributes. So far, I have tried

>> multiple things:

>>

>> 1. The function from this website: http://rstudio-pubs-static.s3.

>> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem

>> with

>> it is that the function cuts the lines in segments of 10km but it is not

>> possible to get the attributes.

>>

>> 2. I have also tried the function from this website:

>> http://stackoverflow.com/questions/38700246/how-do-i-split-

>> divide-polyline-

>> shapefiles-into-equally-length-smaller-segments. The problem is that do

>> not

>> cut the segments in 10km, and it seems to cut the segments with random

>> lengths.

>>

>> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

>> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

>> a7583629 propose to divide lines by segments but do not specify how it is

>> possible to do in based on the length of the segment.

>>

>> Is there any function to do this in R? Or does somebody know the solution

>> for that?

>>

>> My data looks like this:

>>

>> #Script:

>>

>> path_data <- "C:/R/"

>>

>> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>>

>>

>> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklp

>> d3ZnRzA?usp=sharing

>>

>> [[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

### Re: split SLDF with attributes

Marta,

Option 1) using CreateSegments() seems to do what you want. Why don't

you simply use rgeos::over() or rgeos::gTouches() between your old and

new SLDFs to match the index positions of your old attributes to your

new segments?

--Mel.

On 2/23/2017 2:27 PM, marta azores wrote:

> Hi all,

>

> I am working with a large datasets of boat positions, and based on this

> points I created a SpatialLinesDataFrame. Now, I would like to

> split/divide SpatialLinesDataframe

> into n segments of 10km and keep the attributes. So far, I have tried

> multiple things:

>

> 1. The function from this website: http://rstudio-pubs-static.s3.

> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem with

> it is that the function cuts the lines in segments of 10km but it is not

> possible to get the attributes.

>

> 2. I have also tried the function from this website:

> http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

> shapefiles-into-equally-length-smaller-segments. The problem is that do not

> cut the segments in 10km, and it seems to cut the segments with random

> lengths.

>

> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

> a7583629 propose to divide lines by segments but do not specify how it is

> possible to do in based on the length of the segment.

>

> Is there any function to do this in R? Or does somebody know the solution

> for that?

>

> My data looks like this:

>

> #Script:

>

> path_data <- "C:/R/"

>

> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>

>

> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Option 1) using CreateSegments() seems to do what you want. Why don't

you simply use rgeos::over() or rgeos::gTouches() between your old and

new SLDFs to match the index positions of your old attributes to your

new segments?

--Mel.

On 2/23/2017 2:27 PM, marta azores wrote:

> Hi all,

>

> I am working with a large datasets of boat positions, and based on this

> points I created a SpatialLinesDataFrame. Now, I would like to

> split/divide SpatialLinesDataframe

> into n segments of 10km and keep the attributes. So far, I have tried

> multiple things:

>

> 1. The function from this website: http://rstudio-pubs-static.s3.

> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem with

> it is that the function cuts the lines in segments of 10km but it is not

> possible to get the attributes.

>

> 2. I have also tried the function from this website:

> http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

> shapefiles-into-equally-length-smaller-segments. The problem is that do not

> cut the segments in 10km, and it seems to cut the segments with random

> lengths.

>

> 3. On the follwoing website: http://r-sig-geo.2731867.n2.

> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

> a7583629 propose to divide lines by segments but do not specify how it is

> possible to do in based on the length of the segment.

>

> Is there any function to do this in R? Or does somebody know the solution

> for that?

>

> My data looks like this:

>

> #Script:

>

> path_data <- "C:/R/"

>

> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

>

>

> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: calculate selected area from stacked raster object

Thanks for providing an example, see the suggestion inline

On 22/02/2017 10:45, Francesc Montserrat wrote:

> Dear List,

>

> As I have only recently started using R for spatial analysis, and I am

> not a geographer or spatial data specialist by any means, I have a -what

> I presume to be- relatively simple question. I am trying to calculate

> the area of part of a stacked raster object that meets certain

> conditions. More specifically, from a dataset from the deep sea in the

> south Atlantic, I have stacked two raster objects (depth and slope) that

> are further identical in coordinate system (WGS84), x-y (Lat-Long)

> position and projection. From the stacked raster object, I would like to

> extract the part that sits between (say) 1000 and 4000 m depth, with a

> slope of more than 10 degrees. I would like to know what the areal

> extent is in square km and I would like to add it to a previously

> plotted map. Below is a reproducible example:

>

> # Raster object containing depth values

> dpt <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

> ymn=-33.8875, ymx=-28.70417)

> values(dpt) <- sample(-200:-5000, size=(nrow(dpt)*ncol(dpt)), replace=T)

>

> # Raster object containing slope values

> slp <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

> ymn=-33.8875, ymx=-28.70417)

> values(slp) <- sample(0:30, size=(nrow(slp)*ncol(slp)), replace=T)

>

> # Stack raster objects

> stk <- stack(dpt,slp)

> # Define function to be passed to calc

fun <- function(x) {

# consider x to be a vector of length 2

ifelse(x[1] <= -1000 & x[1] >= -4000 & x[2] >= 10, 1, NA)

}

condition_raster <- calc(x = stk, fun = fun)

area_raster <- area(condition_raster, na.rm = TRUE)

cellStats(area_raster, 'sum')

Cheers,

Loïc

>

> # Colour palette

> colrs <-

> colorRampPalette(c("navyblue","dodgerblue3","cyan2","green2","darkgoldenrod1"))

>

> # Plot raster map; does not look like ocean floor because of "sample"

> plot(dpt, xlab="Longitude", ylab="Latitude", col=colrs(100), font.lab=2,

> cex.lab=1.5, las=1)

>

>

> # Create a blank copy of previous raster plot

> selectAtt <- raster(dpt)

> # Fill in cells where Attribute(s) meet(s) conditions

> selectAtt[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >=

> 10] <- 90

> # Set object projection

> projection(selectAtt) <- CRS("+proj=longlat +ellps=WGS84")

> # Plot selection in previous raster

> plot(selectAtt, col="red", add=TRUE, legend=F, proj4string=crswgs84)

>

> #####>>> What is the area (within total area) with both depth (layer.1)

> and slope (layer.2) meeting given conditions ?

>

> a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

> area(a, na.rm=T)

>

> # This gives me the error message:

> Error in (function (classes, fdef, mtable) :

> unable to find an inherited method for function ‘area’ for

> signature ‘"matrix"’

>

> I have tried to find what this actually means, but I cannot find a way

> to solve this. What am I missing here ? Any help is greatly appreciated !

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

On 22/02/2017 10:45, Francesc Montserrat wrote:

> Dear List,

>

> As I have only recently started using R for spatial analysis, and I am

> not a geographer or spatial data specialist by any means, I have a -what

> I presume to be- relatively simple question. I am trying to calculate

> the area of part of a stacked raster object that meets certain

> conditions. More specifically, from a dataset from the deep sea in the

> south Atlantic, I have stacked two raster objects (depth and slope) that

> are further identical in coordinate system (WGS84), x-y (Lat-Long)

> position and projection. From the stacked raster object, I would like to

> extract the part that sits between (say) 1000 and 4000 m depth, with a

> slope of more than 10 degrees. I would like to know what the areal

> extent is in square km and I would like to add it to a previously

> plotted map. Below is a reproducible example:

>

> # Raster object containing depth values

> dpt <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

> ymn=-33.8875, ymx=-28.70417)

> values(dpt) <- sample(-200:-5000, size=(nrow(dpt)*ncol(dpt)), replace=T)

>

> # Raster object containing slope values

> slp <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

> ymn=-33.8875, ymx=-28.70417)

> values(slp) <- sample(0:30, size=(nrow(slp)*ncol(slp)), replace=T)

>

> # Stack raster objects

> stk <- stack(dpt,slp)

> # Define function to be passed to calc

fun <- function(x) {

# consider x to be a vector of length 2

ifelse(x[1] <= -1000 & x[1] >= -4000 & x[2] >= 10, 1, NA)

}

condition_raster <- calc(x = stk, fun = fun)

area_raster <- area(condition_raster, na.rm = TRUE)

cellStats(area_raster, 'sum')

Cheers,

Loïc

>

> # Colour palette

> colrs <-

> colorRampPalette(c("navyblue","dodgerblue3","cyan2","green2","darkgoldenrod1"))

>

> # Plot raster map; does not look like ocean floor because of "sample"

> plot(dpt, xlab="Longitude", ylab="Latitude", col=colrs(100), font.lab=2,

> cex.lab=1.5, las=1)

>

>

> # Create a blank copy of previous raster plot

> selectAtt <- raster(dpt)

> # Fill in cells where Attribute(s) meet(s) conditions

> selectAtt[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >=

> 10] <- 90

> # Set object projection

> projection(selectAtt) <- CRS("+proj=longlat +ellps=WGS84")

> # Plot selection in previous raster

> plot(selectAtt, col="red", add=TRUE, legend=F, proj4string=crswgs84)

>

> #####>>> What is the area (within total area) with both depth (layer.1)

> and slope (layer.2) meeting given conditions ?

>

> a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

> area(a, na.rm=T)

>

> # This gives me the error message:

> Error in (function (classes, fdef, mtable) :

> unable to find an inherited method for function ‘area’ for

> signature ‘"matrix"’

>

> I have tried to find what this actually means, but I cannot find a way

> to solve this. What am I missing here ? Any help is greatly appreciated !

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### split SLDF with attributes

Hi all,

I am working with a large datasets of boat positions, and based on this

points I created a SpatialLinesDataFrame. Now, I would like to

split/divide SpatialLinesDataframe

into n segments of 10km and keep the attributes. So far, I have tried

multiple things:

1. The function from this website: http://rstudio-pubs-static.s3.

amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem with

it is that the function cuts the lines in segments of 10km but it is not

possible to get the attributes.

2. I have also tried the function from this website:

http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

shapefiles-into-equally-length-smaller-segments. The problem is that do not

cut the segments in 10km, and it seems to cut the segments with random

lengths.

3. On the follwoing website: http://r-sig-geo.2731867.n2.

nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

a7583629 propose to divide lines by segments but do not specify how it is

possible to do in based on the length of the segment.

Is there any function to do this in R? Or does somebody know the solution

for that?

My data looks like this:

#Script:

path_data <- "C:/R/"

mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I am working with a large datasets of boat positions, and based on this

points I created a SpatialLinesDataFrame. Now, I would like to

split/divide SpatialLinesDataframe

into n segments of 10km and keep the attributes. So far, I have tried

multiple things:

1. The function from this website: http://rstudio-pubs-static.s3.

amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem with

it is that the function cuts the lines in segments of 10km but it is not

possible to get the attributes.

2. I have also tried the function from this website:

http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-

shapefiles-into-equally-length-smaller-segments. The problem is that do not

cut the segments in 10km, and it seems to cut the segments with random

lengths.

3. On the follwoing website: http://r-sig-geo.2731867.n2.

nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#

a7583629 propose to divide lines by segments but do not specify how it is

possible to do in based on the length of the segment.

Is there any function to do this in R? Or does somebody know the solution

for that?

My data looks like this:

#Script:

path_data <- "C:/R/"

mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")

https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Origin and Destination-specific doubly constrained model

Hello everybody,

I would like to share a difficulty I have that it can be very easy for somebody in this forum.

I have calculated a distance parameter for the whole system in a doubly constrained model :

Model01 <- glm(Flow~Origin+Destination+log(Dij), family=poisson(link="log"),data=Base)

Now I would like to do it, now origin-specific in the sense of Fotheringham (1981, 1983 and so one), that is, for each origin, a distance decay paramater with the purpose to see the spatial structure of the flows.

I suppose it is as bellow meanwhhile I am not sure.

Model02 <- glm(Flow~Origin+Destination+factor(Origin)xlog(Dij), family=poisson(link="log"),data=Base)

Thanks in advance,

Ticiana

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I would like to share a difficulty I have that it can be very easy for somebody in this forum.

I have calculated a distance parameter for the whole system in a doubly constrained model :

Model01 <- glm(Flow~Origin+Destination+log(Dij), family=poisson(link="log"),data=Base)

Now I would like to do it, now origin-specific in the sense of Fotheringham (1981, 1983 and so one), that is, for each origin, a distance decay paramater with the purpose to see the spatial structure of the flows.

I suppose it is as bellow meanwhhile I am not sure.

Model02 <- glm(Flow~Origin+Destination+factor(Origin)xlog(Dij), family=poisson(link="log"),data=Base)

Thanks in advance,

Ticiana

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### calculate selected area from stacked raster object

Dear List,

As I have only recently started using R for spatial analysis, and I am

not a geographer or spatial data specialist by any means, I have a -what

I presume to be- relatively simple question. I am trying to calculate

the area of part of a stacked raster object that meets certain

conditions. More specifically, from a dataset from the deep sea in the

south Atlantic, I have stacked two raster objects (depth and slope) that

are further identical in coordinate system (WGS84), x-y (Lat-Long)

position and projection. From the stacked raster object, I would like to

extract the part that sits between (say) 1000 and 4000 m depth, with a

slope of more than 10 degrees. I would like to know what the areal

extent is in square km and I would like to add it to a previously

plotted map. Below is a reproducible example:

# Raster object containing depth values

dpt <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

ymn=-33.8875, ymx=-28.70417)

values(dpt) <- sample(-200:-5000, size=(nrow(dpt)*ncol(dpt)), replace=T)

# Raster object containing slope values

slp <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

ymn=-33.8875, ymx=-28.70417)

values(slp) <- sample(0:30, size=(nrow(slp)*ncol(slp)), replace=T)

# Stack raster objects

stk <- stack(dpt,slp)

# Colour palette

colrs <-

colorRampPalette(c("navyblue","dodgerblue3","cyan2","green2","darkgoldenrod1"))

# Plot raster map; does not look like ocean floor because of "sample"

plot(dpt, xlab="Longitude", ylab="Latitude", col=colrs(100), font.lab=2,

cex.lab=1.5, las=1)

# Create a blank copy of previous raster plot

selectAtt <- raster(dpt)

# Fill in cells where Attribute(s) meet(s) conditions

selectAtt[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >=

10] <- 90

# Set object projection

projection(selectAtt) <- CRS("+proj=longlat +ellps=WGS84")

# Plot selection in previous raster

plot(selectAtt, col="red", add=TRUE, legend=F, proj4string=crswgs84)

#####>>> What is the area (within total area) with both depth (layer.1)

and slope (layer.2) meeting given conditions ?

a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

area(a, na.rm=T)

# This gives me the error message:

Error in (function (classes, fdef, mtable) :

unable to find an inherited method for function ‘area’ for

signature ‘"matrix"’

I have tried to find what this actually means, but I cannot find a way

to solve this. What am I missing here ? Any help is greatly appreciated !

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

As I have only recently started using R for spatial analysis, and I am

not a geographer or spatial data specialist by any means, I have a -what

I presume to be- relatively simple question. I am trying to calculate

the area of part of a stacked raster object that meets certain

conditions. More specifically, from a dataset from the deep sea in the

south Atlantic, I have stacked two raster objects (depth and slope) that

are further identical in coordinate system (WGS84), x-y (Lat-Long)

position and projection. From the stacked raster object, I would like to

extract the part that sits between (say) 1000 and 4000 m depth, with a

slope of more than 10 degrees. I would like to know what the areal

extent is in square km and I would like to add it to a previously

plotted map. Below is a reproducible example:

# Raster object containing depth values

dpt <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

ymn=-33.8875, ymx=-28.70417)

values(dpt) <- sample(-200:-5000, size=(nrow(dpt)*ncol(dpt)), replace=T)

# Raster object containing slope values

slp <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,

ymn=-33.8875, ymx=-28.70417)

values(slp) <- sample(0:30, size=(nrow(slp)*ncol(slp)), replace=T)

# Stack raster objects

stk <- stack(dpt,slp)

# Colour palette

colrs <-

colorRampPalette(c("navyblue","dodgerblue3","cyan2","green2","darkgoldenrod1"))

# Plot raster map; does not look like ocean floor because of "sample"

plot(dpt, xlab="Longitude", ylab="Latitude", col=colrs(100), font.lab=2,

cex.lab=1.5, las=1)

# Create a blank copy of previous raster plot

selectAtt <- raster(dpt)

# Fill in cells where Attribute(s) meet(s) conditions

selectAtt[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >=

10] <- 90

# Set object projection

projection(selectAtt) <- CRS("+proj=longlat +ellps=WGS84")

# Plot selection in previous raster

plot(selectAtt, col="red", add=TRUE, legend=F, proj4string=crswgs84)

#####>>> What is the area (within total area) with both depth (layer.1)

and slope (layer.2) meeting given conditions ?

a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

area(a, na.rm=T)

# This gives me the error message:

Error in (function (classes, fdef, mtable) :

unable to find an inherited method for function ‘area’ for

signature ‘"matrix"’

I have tried to find what this actually means, but I cannot find a way

to solve this. What am I missing here ? Any help is greatly appreciated !

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: Subset a SpLineDF

By assigning slots directly, you always take a serious risk that you

create objects that are considered invalid in a next step, as is the

case in your example.

My guess is that when you subset @data, you should also replace the

@lines slot with the corresponding subset of @lines (which is a list).

The more straightforward (and safe) way is to subset the whole object as

if it were a data.frame:

TRONCON_TOPO_DSAV_coor[! TRONCON_TOPO_DSAV_coord$doublon, ]

(note the absence of @)

On 22/02/17 11:49, Tristan Bourgeois wrote:

> Dear all,

>

> I'm facing a problem I can't solve even if the solution should not be so

> hard to find.

>

> I want to delete the duplicates in a spatial line dataframe .Those

> duplicates have been created by shapesfiles merging on QGis (boundary

> rivers are duplicated by the merging of neighbouring departments layers)

>

> So here's the code I wrote :

>

>

> TRONCON_TOPO_DSAV_coord <-readOGR(dsn = "Sorties/Tronçon

> BDTOPO/TRONCON_TOPO_DSAV/Avec coordonnées", layer =

> "TRONCON_TOPO_DSAV_coord")

> TRONCON_TOPO_DSAV_coord@data

> $doublon<-duplicated(TRONCON_TOPO_DSAV_coord@data$ID)### Create a fied for

> duplicate identification

> TRONCON_TOPO_DSAV_coord@data<-subset(TRONCON_TOPO_DSAV_coord@data

> ,TRONCON_TOPO_DSAV_coord@data$doublon=="FALSE")### Subset to delete the

> duplicate

> TRONCON_TOPO_DSAV_coord@data ### Checking results

>

> writeOGR(TRONCON_TOPO_DSAV_coord,dsn =

> path,layer="TRONCON_TOPO_DSAV_coord_unique_test",driver ="ESRI Shapefile")

>

>

> When I want to export my OGR file, I get this error message :"Error in

> writeOGR(TRONCON_TOPO_DSAV_coord, dsn = path, layer =

> "TRONCON_TOPO_DSAV_coord_unique_test", :

> number of objects mismatch

>

>

> It seems that I have not the same number of objects in the different

> shapefile slots. How to update these objets after my subset ?

>

>

> Cheers !

>

> "

> --

Edzer Pebesma

Institute for Geoinformatics (ifgi), University of Münster

Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081

Journal of Statistical Software: http://www.jstatsoft.org/

Computers & Geosciences: http://elsevier.com/locate/cageo/

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

create objects that are considered invalid in a next step, as is the

case in your example.

My guess is that when you subset @data, you should also replace the

@lines slot with the corresponding subset of @lines (which is a list).

The more straightforward (and safe) way is to subset the whole object as

if it were a data.frame:

TRONCON_TOPO_DSAV_coor[! TRONCON_TOPO_DSAV_coord$doublon, ]

(note the absence of @)

On 22/02/17 11:49, Tristan Bourgeois wrote:

> Dear all,

>

> I'm facing a problem I can't solve even if the solution should not be so

> hard to find.

>

> I want to delete the duplicates in a spatial line dataframe .Those

> duplicates have been created by shapesfiles merging on QGis (boundary

> rivers are duplicated by the merging of neighbouring departments layers)

>

> So here's the code I wrote :

>

>

> TRONCON_TOPO_DSAV_coord <-readOGR(dsn = "Sorties/Tronçon

> BDTOPO/TRONCON_TOPO_DSAV/Avec coordonnées", layer =

> "TRONCON_TOPO_DSAV_coord")

> TRONCON_TOPO_DSAV_coord@data

> $doublon<-duplicated(TRONCON_TOPO_DSAV_coord@data$ID)### Create a fied for

> duplicate identification

> TRONCON_TOPO_DSAV_coord@data<-subset(TRONCON_TOPO_DSAV_coord@data

> ,TRONCON_TOPO_DSAV_coord@data$doublon=="FALSE")### Subset to delete the

> duplicate

> TRONCON_TOPO_DSAV_coord@data ### Checking results

>

> writeOGR(TRONCON_TOPO_DSAV_coord,dsn =

> path,layer="TRONCON_TOPO_DSAV_coord_unique_test",driver ="ESRI Shapefile")

>

>

> When I want to export my OGR file, I get this error message :"Error in

> writeOGR(TRONCON_TOPO_DSAV_coord, dsn = path, layer =

> "TRONCON_TOPO_DSAV_coord_unique_test", :

> number of objects mismatch

>

>

> It seems that I have not the same number of objects in the different

> shapefile slots. How to update these objets after my subset ?

>

>

> Cheers !

>

> "

> --

Edzer Pebesma

Institute for Geoinformatics (ifgi), University of Münster

Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081

Journal of Statistical Software: http://www.jstatsoft.org/

Computers & Geosciences: http://elsevier.com/locate/cageo/

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

**signature.asc**(484 bytes) Download Attachment### Subset a SpLineDF

Dear all,

I'm facing a problem I can't solve even if the solution should not be so

hard to find.

I want to delete the duplicates in a spatial line dataframe .Those

duplicates have been created by shapesfiles merging on QGis (boundary

rivers are duplicated by the merging of neighbouring departments layers)

So here's the code I wrote :

TRONCON_TOPO_DSAV_coord <-readOGR(dsn = "Sorties/Tronçon

BDTOPO/TRONCON_TOPO_DSAV/Avec coordonnées", layer =

"TRONCON_TOPO_DSAV_coord")

TRONCON_TOPO_DSAV_coord@data

$doublon<-duplicated(TRONCON_TOPO_DSAV_coord@data$ID)### Create a fied for

duplicate identification

TRONCON_TOPO_DSAV_coord@data<-subset(TRONCON_TOPO_DSAV_coord@data

,TRONCON_TOPO_DSAV_coord@data$doublon=="FALSE")### Subset to delete the

duplicate

TRONCON_TOPO_DSAV_coord@data ### Checking results

writeOGR(TRONCON_TOPO_DSAV_coord,dsn =

path,layer="TRONCON_TOPO_DSAV_coord_unique_test",driver ="ESRI Shapefile")

When I want to export my OGR file, I get this error message :"Error in

writeOGR(TRONCON_TOPO_DSAV_coord, dsn = path, layer =

"TRONCON_TOPO_DSAV_coord_unique_test", :

number of objects mismatch

It seems that I have not the same number of objects in the different

shapefile slots. How to update these objets after my subset ?

Cheers !

"

--

Tristan Bourgeois

Hydromorphologist

Seine Normandie Water Agency

Rouen, France

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I'm facing a problem I can't solve even if the solution should not be so

hard to find.

I want to delete the duplicates in a spatial line dataframe .Those

duplicates have been created by shapesfiles merging on QGis (boundary

rivers are duplicated by the merging of neighbouring departments layers)

So here's the code I wrote :

TRONCON_TOPO_DSAV_coord <-readOGR(dsn = "Sorties/Tronçon

BDTOPO/TRONCON_TOPO_DSAV/Avec coordonnées", layer =

"TRONCON_TOPO_DSAV_coord")

TRONCON_TOPO_DSAV_coord@data

$doublon<-duplicated(TRONCON_TOPO_DSAV_coord@data$ID)### Create a fied for

duplicate identification

TRONCON_TOPO_DSAV_coord@data<-subset(TRONCON_TOPO_DSAV_coord@data

,TRONCON_TOPO_DSAV_coord@data$doublon=="FALSE")### Subset to delete the

duplicate

TRONCON_TOPO_DSAV_coord@data ### Checking results

writeOGR(TRONCON_TOPO_DSAV_coord,dsn =

path,layer="TRONCON_TOPO_DSAV_coord_unique_test",driver ="ESRI Shapefile")

When I want to export my OGR file, I get this error message :"Error in

writeOGR(TRONCON_TOPO_DSAV_coord, dsn = path, layer =

"TRONCON_TOPO_DSAV_coord_unique_test", :

number of objects mismatch

It seems that I have not the same number of objects in the different

shapefile slots. How to update these objets after my subset ?

Cheers !

"

--

Tristan Bourgeois

Hydromorphologist

Seine Normandie Water Agency

Rouen, France

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: Optimizing spatial query in R

I generally use postgis for this sort of thing.

THK

http://www.keittlab.org/

On Mon, Feb 20, 2017 at 2:45 AM, Ervan Rutishauser <[hidden email]

> wrote:

> Dear All,

>

> I am desperately trying to fasten my algorithm to estimate the fraction of

> tree crown that overlap a given 10x10 subplot in a forest plot. I have

> combined a set of spatial functions (gDistance, extract) & objects

> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

> forest plot).

>

> My detailed problem and a reproducible example are posted on Stackoverflow

> <http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r>

> and

> append below, if you wanna have a look. Apart of Amazon Web Server, is

> anyone aware of a sever where to execute (and save results back) R codes

> online?

> Thanks for any help.

> Best regards,

>

>

> # A. Define objects

> require(sp)

> require(raster)

> require(rgdal)

> require(rgeos)

> require(dismo)

> radius=25 # max search radius around 10 x 10 m cells

> res <- vector() # where to store results

>

> # Create a fake set of trees with x,y coordinates and trunk diameter

> (=dbh)

> set.seed(0)

> survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,

> replace=T),dbh=sample(100,1000,replace=T))

> coordinates(survey) <- ~x+y

>

> # Define 10 x 10 subplots

> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

> survey$subplot <- over(survey,grid10)

> # B. Now find fraction of tree crown overlapping each subplot

> for (i in 1:100) {

> # Extract centroïd of each the ith cell

> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

> corner <-

> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$

> x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>

>

> # Find trees in a max radius (define above)

> tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=

> radius^2),]

>

>

> # Define tree crown based on tree diameter

> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

> meter

>

> # Compute the distance from each tree to cell's borders

> pDist <- vector()

> for (k in 1:nrow(tem)) {

> pDist[k] <-

> gDistance(tem[k,],SpatialPolygons(list(Polygons(

> list(Polygon(corner)),1))))

> }

>

> # Keeps only trees whose crown is lower than the above

> distance (=overlap)

> overlap.trees <- tem[which(pDist<=tem$crownr),]

> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown

> area

>

> # Creat polygons from overlapping crowns

> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

> lonlat=F, dissolve=F)

> crown <- polygons(c1)

> Crown <-

> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(

> dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

>

> # Create a fine grid points to retrieve the fraction of

> overlapping crowns

> max.dist <- ceiling(sqrt(which.max((centro$x -

> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

> to narrow search

>

> finegrid <-

> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$

> x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

> coordinates(finegrid) <- ~ x+y

> A <- extract(Crown,finegrid)

> Crown@data$ID <- seq(1,length(crown),1)

> B <- as.data.frame(table(A$poly.ID))

> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

> B$overlap <- B$Freq/B$crown.area

> B$overlap[B$overlap>1] <- 1

> res[i] <- sum(B$overlap)

> }

> # C. Check the result

> res # sum of crown fraction overlapping each cell (works fine)

>

> --

> _____________________________

> Ervan Rutishauser, PhD

> <http://carboforexpert.ch/>STRI post-doctoral fellow

> CarboForExpert

> <http://carboforexpert.ch/>

> Skype: ervan-CH

>

> [[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

THK

http://www.keittlab.org/

On Mon, Feb 20, 2017 at 2:45 AM, Ervan Rutishauser <[hidden email]

> wrote:

> Dear All,

>

> I am desperately trying to fasten my algorithm to estimate the fraction of

> tree crown that overlap a given 10x10 subplot in a forest plot. I have

> combined a set of spatial functions (gDistance, extract) & objects

> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

> forest plot).

>

> My detailed problem and a reproducible example are posted on Stackoverflow

> <http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r>

> and

> append below, if you wanna have a look. Apart of Amazon Web Server, is

> anyone aware of a sever where to execute (and save results back) R codes

> online?

> Thanks for any help.

> Best regards,

>

>

> # A. Define objects

> require(sp)

> require(raster)

> require(rgdal)

> require(rgeos)

> require(dismo)

> radius=25 # max search radius around 10 x 10 m cells

> res <- vector() # where to store results

>

> # Create a fake set of trees with x,y coordinates and trunk diameter

> (=dbh)

> set.seed(0)

> survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,

> replace=T),dbh=sample(100,1000,replace=T))

> coordinates(survey) <- ~x+y

>

> # Define 10 x 10 subplots

> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

> survey$subplot <- over(survey,grid10)

> # B. Now find fraction of tree crown overlapping each subplot

> for (i in 1:100) {

> # Extract centroïd of each the ith cell

> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

> corner <-

> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$

> x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>

>

> # Find trees in a max radius (define above)

> tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=

> radius^2),]

>

>

> # Define tree crown based on tree diameter

> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

> meter

>

> # Compute the distance from each tree to cell's borders

> pDist <- vector()

> for (k in 1:nrow(tem)) {

> pDist[k] <-

> gDistance(tem[k,],SpatialPolygons(list(Polygons(

> list(Polygon(corner)),1))))

> }

>

> # Keeps only trees whose crown is lower than the above

> distance (=overlap)

> overlap.trees <- tem[which(pDist<=tem$crownr),]

> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown

> area

>

> # Creat polygons from overlapping crowns

> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

> lonlat=F, dissolve=F)

> crown <- polygons(c1)

> Crown <-

> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(

> dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

>

> # Create a fine grid points to retrieve the fraction of

> overlapping crowns

> max.dist <- ceiling(sqrt(which.max((centro$x -

> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

> to narrow search

>

> finegrid <-

> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$

> x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

> coordinates(finegrid) <- ~ x+y

> A <- extract(Crown,finegrid)

> Crown@data$ID <- seq(1,length(crown),1)

> B <- as.data.frame(table(A$poly.ID))

> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

> B$overlap <- B$Freq/B$crown.area

> B$overlap[B$overlap>1] <- 1

> res[i] <- sum(B$overlap)

> }

> # C. Check the result

> res # sum of crown fraction overlapping each cell (works fine)

>

> --

> _____________________________

> Ervan Rutishauser, PhD

> <http://carboforexpert.ch/>STRI post-doctoral fellow

> CarboForExpert

> <http://carboforexpert.ch/>

> Skype: ervan-CH

>

> [[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

### Re: Optimizing spatial query in R

The biggest bottleneck was creating a SpatialPolygons object and

calculating gDistance. After computing the SP object outside the for loop,

the run time came down dramatically (see SO post).

Cheers,

Roman

On Mon, Feb 20, 2017 at 10:56 AM, Chris Reudenbach <

[hidden email]> wrote:

> Ervan,

>

> Just on the run, especially if you deal with real data and some 100K

> -1000K of vectors I think the fastest way to do so is to engage GDAL

> GRASS7/SAGA and their ability to deal highly efficient with this kind of

> topological and geometrical queries.

>

> A very good pure R alternative is to use the sf package that is ways

> faster and more satble than sp. It also provides this type of GIS

> functionality like st_overlaps() etc.

>

> cheers Chris

>

>

>

>

> On 20.02.2017 09:45, Ervan Rutishauser wrote:

>

>> Dear All,

>>

>> I am desperately trying to fasten my algorithm to estimate the fraction of

>> tree crown that overlap a given 10x10 subplot in a forest plot. I have

>> combined a set of spatial functions (gDistance, extract) & objects

>> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

>> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

>> forest plot).

>>

>> My detailed problem and a reproducible example are posted on Stackoverflow

>> <http://stackoverflow.com/questions/42303559/optimizing-spat

>> ial-query-in-r> and

>> append below, if you wanna have a look. Apart of Amazon Web Server, is

>> anyone aware of a sever where to execute (and save results back) R codes

>> online?

>> Thanks for any help.

>> Best regards,

>>

>>

>> # A. Define objects

>> require(sp)

>> require(raster)

>> require(rgdal)

>> require(rgeos)

>> require(dismo)

>> radius=25 # max search radius around 10 x 10 m cells

>> res <- vector() # where to store results

>>

>> # Create a fake set of trees with x,y coordinates and trunk diameter

>> (=dbh)

>> set.seed(0)

>> survey <- data.frame(x=sample(99,1000,re

>> place=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

>> coordinates(survey) <- ~x+y

>>

>> # Define 10 x 10 subplots

>> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

>> survey$subplot <- over(survey,grid10)

>> # B. Now find fraction of tree crown overlapping each subplot

>> for (i in 1:100) {

>> # Extract centroïd of each the ith cell

>> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

>> corner <-

>> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),

>> y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>>

>>

>> # Find trees in a max radius (define above)

>> tem <- survey[which((centro$x-survey$

>> x)^2+(centro$y-survey$y)^2<=radius^2),]

>>

>>

>> # Define tree crown based on tree diameter

>> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

>> meter

>>

>> # Compute the distance from each tree to cell's borders

>> pDist <- vector()

>> for (k in 1:nrow(tem)) {

>> pDist[k] <-

>> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(

>> Polygon(corner)),1))))

>> }

>>

>> # Keeps only trees whose crown is lower than the above

>> distance (=overlap)

>> overlap.trees <- tem[which(pDist<=tem$crownr),]

>> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute

>> crown area

>>

>> # Creat polygons from overlapping crowns

>> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

>> lonlat=F, dissolve=F)

>> crown <- polygons(c1)

>> Crown <-

>> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=

>> overlap.trees$dbh,crown.area=overlap.trees$crowna))

>>

>> # Create a fine grid points to retrieve the fraction of

>> overlapping crowns

>> max.dist <- ceiling(sqrt(which.max((centro$x -

>> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

>> to narrow search

>>

>> finegrid <-

>> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+

>> max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

>> coordinates(finegrid) <- ~ x+y

>> A <- extract(Crown,finegrid)

>> Crown@data$ID <- seq(1,length(crown),1)

>> B <- as.data.frame(table(A$poly.ID))

>> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

>> B$overlap <- B$Freq/B$crown.area

>> B$overlap[B$overlap>1] <- 1

>> res[i] <- sum(B$overlap)

>> }

>> # C. Check the result

>> res # sum of crown fraction overlapping each cell (works fine)

>>

>>

> --

> Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

> Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

> Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

> gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

>

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

--

In God we trust, all others bring data.

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

calculating gDistance. After computing the SP object outside the for loop,

the run time came down dramatically (see SO post).

Cheers,

Roman

On Mon, Feb 20, 2017 at 10:56 AM, Chris Reudenbach <

[hidden email]> wrote:

> Ervan,

>

> Just on the run, especially if you deal with real data and some 100K

> -1000K of vectors I think the fastest way to do so is to engage GDAL

> GRASS7/SAGA and their ability to deal highly efficient with this kind of

> topological and geometrical queries.

>

> A very good pure R alternative is to use the sf package that is ways

> faster and more satble than sp. It also provides this type of GIS

> functionality like st_overlaps() etc.

>

> cheers Chris

>

>

>

>

> On 20.02.2017 09:45, Ervan Rutishauser wrote:

>

>> Dear All,

>>

>> I am desperately trying to fasten my algorithm to estimate the fraction of

>> tree crown that overlap a given 10x10 subplot in a forest plot. I have

>> combined a set of spatial functions (gDistance, extract) & objects

>> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

>> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

>> forest plot).

>>

>> My detailed problem and a reproducible example are posted on Stackoverflow

>> <http://stackoverflow.com/questions/42303559/optimizing-spat

>> ial-query-in-r> and

>> append below, if you wanna have a look. Apart of Amazon Web Server, is

>> anyone aware of a sever where to execute (and save results back) R codes

>> online?

>> Thanks for any help.

>> Best regards,

>>

>>

>> # A. Define objects

>> require(sp)

>> require(raster)

>> require(rgdal)

>> require(rgeos)

>> require(dismo)

>> radius=25 # max search radius around 10 x 10 m cells

>> res <- vector() # where to store results

>>

>> # Create a fake set of trees with x,y coordinates and trunk diameter

>> (=dbh)

>> set.seed(0)

>> survey <- data.frame(x=sample(99,1000,re

>> place=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

>> coordinates(survey) <- ~x+y

>>

>> # Define 10 x 10 subplots

>> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

>> survey$subplot <- over(survey,grid10)

>> # B. Now find fraction of tree crown overlapping each subplot

>> for (i in 1:100) {

>> # Extract centroïd of each the ith cell

>> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

>> corner <-

>> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),

>> y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>>

>>

>> # Find trees in a max radius (define above)

>> tem <- survey[which((centro$x-survey$

>> x)^2+(centro$y-survey$y)^2<=radius^2),]

>>

>>

>> # Define tree crown based on tree diameter

>> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

>> meter

>>

>> # Compute the distance from each tree to cell's borders

>> pDist <- vector()

>> for (k in 1:nrow(tem)) {

>> pDist[k] <-

>> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(

>> Polygon(corner)),1))))

>> }

>>

>> # Keeps only trees whose crown is lower than the above

>> distance (=overlap)

>> overlap.trees <- tem[which(pDist<=tem$crownr),]

>> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute

>> crown area

>>

>> # Creat polygons from overlapping crowns

>> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

>> lonlat=F, dissolve=F)

>> crown <- polygons(c1)

>> Crown <-

>> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=

>> overlap.trees$dbh,crown.area=overlap.trees$crowna))

>>

>> # Create a fine grid points to retrieve the fraction of

>> overlapping crowns

>> max.dist <- ceiling(sqrt(which.max((centro$x -

>> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

>> to narrow search

>>

>> finegrid <-

>> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+

>> max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

>> coordinates(finegrid) <- ~ x+y

>> A <- extract(Crown,finegrid)

>> Crown@data$ID <- seq(1,length(crown),1)

>> B <- as.data.frame(table(A$poly.ID))

>> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

>> B$overlap <- B$Freq/B$crown.area

>> B$overlap[B$overlap>1] <- 1

>> res[i] <- sum(B$overlap)

>> }

>> # C. Check the result

>> res # sum of crown fraction overlapping each cell (works fine)

>>

>>

> --

> Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

> Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

> Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

> gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

>

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

--

In God we trust, all others bring data.

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: Optimizing spatial query in R

Dear Chris,

Thanks for your prompt reply and pointing out the "sf" package . I wasn't

aware of it, but I am glad to see it has been released only 2 weeks ago ;)

I'll have a look right now.

Tschüss

Ervan

On 20 February 2017 at 10:56, Chris Reudenbach <[hidden email]>

wrote:

> Ervan,

>

> Just on the run, especially if you deal with real data and some 100K

> -1000K of vectors I think the fastest way to do so is to engage GDAL

> GRASS7/SAGA and their ability to deal highly efficient with this kind of

> topological and geometrical queries.

>

> A very good pure R alternative is to use the sf package that is ways

> faster and more satble than sp. It also provides this type of GIS

> functionality like st_overlaps() etc.

>

> cheers Chris

>

>

>

> On 20.02.2017 09:45, Ervan Rutishauser wrote:

>

>> Dear All,

>>

>> I am desperately trying to fasten my algorithm to estimate the fraction of

>> tree crown that overlap a given 10x10 subplot in a forest plot. I have

>> combined a set of spatial functions (gDistance, extract) & objects

>> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

>> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

>> forest plot).

>>

>> My detailed problem and a reproducible example are posted on Stackoverflow

>> <http://stackoverflow.com/questions/42303559/optimizing-spat

>> ial-query-in-r> and

>>

>> append below, if you wanna have a look. Apart of Amazon Web Server, is

>> anyone aware of a sever where to execute (and save results back) R codes

>> online?

>> Thanks for any help.

>> Best regards,

>>

>>

>> # A. Define objects

>> require(sp)

>> require(raster)

>> require(rgdal)

>> require(rgeos)

>> require(dismo)

>> radius=25 # max search radius around 10 x 10 m cells

>> res <- vector() # where to store results

>>

>> # Create a fake set of trees with x,y coordinates and trunk diameter

>> (=dbh)

>> set.seed(0)

>> survey <- data.frame(x=sample(99,1000,re

>> place=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

>> coordinates(survey) <- ~x+y

>>

>> # Define 10 x 10 subplots

>> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

>> survey$subplot <- over(survey,grid10)

>> # B. Now find fraction of tree crown overlapping each subplot

>> for (i in 1:100) {

>> # Extract centroïd of each the ith cell

>> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

>> corner <-

>> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),

>> y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>>

>>

>> # Find trees in a max radius (define above)

>> tem <- survey[which((centro$x-survey$

>> x)^2+(centro$y-survey$y)^2<=radius^2),]

>>

>>

>> # Define tree crown based on tree diameter

>> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

>> meter

>>

>> # Compute the distance from each tree to cell's borders

>> pDist <- vector()

>> for (k in 1:nrow(tem)) {

>> pDist[k] <-

>> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(

>> Polygon(corner)),1))))

>> }

>>

>> # Keeps only trees whose crown is lower than the above

>> distance (=overlap)

>> overlap.trees <- tem[which(pDist<=tem$crownr),]

>> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute

>> crown area

>>

>> # Creat polygons from overlapping crowns

>> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

>> lonlat=F, dissolve=F)

>> crown <- polygons(c1)

>> Crown <-

>> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=

>> overlap.trees$dbh,crown.area=overlap.trees$crowna))

>>

>> # Create a fine grid points to retrieve the fraction of

>> overlapping crowns

>> max.dist <- ceiling(sqrt(which.max((centro$x -

>> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

>> to narrow search

>>

>> finegrid <-

>> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+

>> max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

>> coordinates(finegrid) <- ~ x+y

>> A <- extract(Crown,finegrid)

>> Crown@data$ID <- seq(1,length(crown),1)

>> B <- as.data.frame(table(A$poly.ID))

>> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

>> B$overlap <- B$Freq/B$crown.area

>> B$overlap[B$overlap>1] <- 1

>> res[i] <- sum(B$overlap)

>> }

>> # C. Check the result

>> res # sum of crown fraction overlapping each cell (works fine)

>>

>>

> --

> Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

> Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

> Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

> gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

>

>

--

_____________________________

Ervan Rutishauser, PhD

<http://carboforexpert.ch/>STRI post-doctoral fellow

CarboForExpert

<http://carboforexpert.ch/>

Skype: ervan-CH

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Thanks for your prompt reply and pointing out the "sf" package . I wasn't

aware of it, but I am glad to see it has been released only 2 weeks ago ;)

I'll have a look right now.

Tschüss

Ervan

On 20 February 2017 at 10:56, Chris Reudenbach <[hidden email]>

wrote:

> Ervan,

>

> Just on the run, especially if you deal with real data and some 100K

> -1000K of vectors I think the fastest way to do so is to engage GDAL

> GRASS7/SAGA and their ability to deal highly efficient with this kind of

> topological and geometrical queries.

>

> A very good pure R alternative is to use the sf package that is ways

> faster and more satble than sp. It also provides this type of GIS

> functionality like st_overlaps() etc.

>

> cheers Chris

>

>

>

> On 20.02.2017 09:45, Ervan Rutishauser wrote:

>

>> Dear All,

>>

>> I am desperately trying to fasten my algorithm to estimate the fraction of

>> tree crown that overlap a given 10x10 subplot in a forest plot. I have

>> combined a set of spatial functions (gDistance, extract) & objects

>> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

>> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

>> forest plot).

>>

>> My detailed problem and a reproducible example are posted on Stackoverflow

>> <http://stackoverflow.com/questions/42303559/optimizing-spat

>> ial-query-in-r> and

>>

>> append below, if you wanna have a look. Apart of Amazon Web Server, is

>> anyone aware of a sever where to execute (and save results back) R codes

>> online?

>> Thanks for any help.

>> Best regards,

>>

>>

>> # A. Define objects

>> require(sp)

>> require(raster)

>> require(rgdal)

>> require(rgeos)

>> require(dismo)

>> radius=25 # max search radius around 10 x 10 m cells

>> res <- vector() # where to store results

>>

>> # Create a fake set of trees with x,y coordinates and trunk diameter

>> (=dbh)

>> set.seed(0)

>> survey <- data.frame(x=sample(99,1000,re

>> place=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

>> coordinates(survey) <- ~x+y

>>

>> # Define 10 x 10 subplots

>> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

>> survey$subplot <- over(survey,grid10)

>> # B. Now find fraction of tree crown overlapping each subplot

>> for (i in 1:100) {

>> # Extract centroïd of each the ith cell

>> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

>> corner <-

>> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),

>> y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>>

>>

>> # Find trees in a max radius (define above)

>> tem <- survey[which((centro$x-survey$

>> x)^2+(centro$y-survey$y)^2<=radius^2),]

>>

>>

>> # Define tree crown based on tree diameter

>> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in

>> meter

>>

>> # Compute the distance from each tree to cell's borders

>> pDist <- vector()

>> for (k in 1:nrow(tem)) {

>> pDist[k] <-

>> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(

>> Polygon(corner)),1))))

>> }

>>

>> # Keeps only trees whose crown is lower than the above

>> distance (=overlap)

>> overlap.trees <- tem[which(pDist<=tem$crownr),]

>> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute

>> crown area

>>

>> # Creat polygons from overlapping crowns

>> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

>> lonlat=F, dissolve=F)

>> crown <- polygons(c1)

>> Crown <-

>> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=

>> overlap.trees$dbh,crown.area=overlap.trees$crowna))

>>

>> # Create a fine grid points to retrieve the fraction of

>> overlapping crowns

>> max.dist <- ceiling(sqrt(which.max((centro$x -

>> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

>> to narrow search

>>

>> finegrid <-

>> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+

>> max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

>> coordinates(finegrid) <- ~ x+y

>> A <- extract(Crown,finegrid)

>> Crown@data$ID <- seq(1,length(crown),1)

>> B <- as.data.frame(table(A$poly.ID))

>> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

>> B$overlap <- B$Freq/B$crown.area

>> B$overlap[B$overlap>1] <- 1

>> res[i] <- sum(B$overlap)

>> }

>> # C. Check the result

>> res # sum of crown fraction overlapping each cell (works fine)

>>

>>

> --

> Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

> Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

> Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

> gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

>

>

--

_____________________________

Ervan Rutishauser, PhD

<http://carboforexpert.ch/>STRI post-doctoral fellow

CarboForExpert

<http://carboforexpert.ch/>

Skype: ervan-CH

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Re: Optimizing spatial query in R

Ervan,

Just on the run, especially if you deal with real data and some 100K

-1000K of vectors I think the fastest way to do so is to engage GDAL

GRASS7/SAGA and their ability to deal highly efficient with this kind

of topological and geometrical queries.

A very good pure R alternative is to use the sf package that is ways

faster and more satble than sp. It also provides this type of GIS

functionality like st_overlaps() etc.

cheers Chris

On 20.02.2017 09:45, Ervan Rutishauser wrote:

> Dear All,

>

> I am desperately trying to fasten my algorithm to estimate the fraction of

> tree crown that overlap a given 10x10 subplot in a forest plot. I have

> combined a set of spatial functions (gDistance, extract) & objects

> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

> forest plot).

>

> My detailed problem and a reproducible example are posted on Stackoverflow

> <http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r> and

> append below, if you wanna have a look. Apart of Amazon Web Server, is

> anyone aware of a sever where to execute (and save results back) R codes

> online?

> Thanks for any help.

> Best regards,

>

>

> # A. Define objects

> require(sp)

> require(raster)

> require(rgdal)

> require(rgeos)

> require(dismo)

> radius=25 # max search radius around 10 x 10 m cells

> res <- vector() # where to store results

>

> # Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)

> set.seed(0)

> survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

> coordinates(survey) <- ~x+y

>

> # Define 10 x 10 subplots

> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

> survey$subplot <- over(survey,grid10)

> # B. Now find fraction of tree crown overlapping each subplot

> for (i in 1:100) {

> # Extract centroïd of each the ith cell

> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

> corner <-

> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>

>

> # Find trees in a max radius (define above)

> tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]

>

>

> # Define tree crown based on tree diameter

> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

>

> # Compute the distance from each tree to cell's borders

> pDist <- vector()

> for (k in 1:nrow(tem)) {

> pDist[k] <-

> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))

> }

>

> # Keeps only trees whose crown is lower than the above

> distance (=overlap)

> overlap.trees <- tem[which(pDist<=tem$crownr),]

> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown area

>

> # Creat polygons from overlapping crowns

> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

> lonlat=F, dissolve=F)

> crown <- polygons(c1)

> Crown <-

> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

>

> # Create a fine grid points to retrieve the fraction of

> overlapping crowns

> max.dist <- ceiling(sqrt(which.max((centro$x -

> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

> to narrow search

>

> finegrid <-

> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

> coordinates(finegrid) <- ~ x+y

> A <- extract(Crown,finegrid)

> Crown@data$ID <- seq(1,length(crown),1)

> B <- as.data.frame(table(A$poly.ID))

> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

> B$overlap <- B$Freq/B$crown.area

> B$overlap[B$overlap>1] <- 1

> res[i] <- sum(B$overlap)

> }

> # C. Check the result

> res # sum of crown fraction overlapping each cell (works fine)

>

--

Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Just on the run, especially if you deal with real data and some 100K

-1000K of vectors I think the fastest way to do so is to engage GDAL

GRASS7/SAGA and their ability to deal highly efficient with this kind

of topological and geometrical queries.

A very good pure R alternative is to use the sf package that is ways

faster and more satble than sp. It also provides this type of GIS

functionality like st_overlaps() etc.

cheers Chris

On 20.02.2017 09:45, Ervan Rutishauser wrote:

> Dear All,

>

> I am desperately trying to fasten my algorithm to estimate the fraction of

> tree crown that overlap a given 10x10 subplot in a forest plot. I have

> combined a set of spatial functions (gDistance, extract) & objects

> (SpatialGrid, SpatialPolygons) in a way that is probably not the most

> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

> forest plot).

>

> My detailed problem and a reproducible example are posted on Stackoverflow

> <http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r> and

> append below, if you wanna have a look. Apart of Amazon Web Server, is

> anyone aware of a sever where to execute (and save results back) R codes

> online?

> Thanks for any help.

> Best regards,

>

>

> # A. Define objects

> require(sp)

> require(raster)

> require(rgdal)

> require(rgeos)

> require(dismo)

> radius=25 # max search radius around 10 x 10 m cells

> res <- vector() # where to store results

>

> # Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)

> set.seed(0)

> survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

> coordinates(survey) <- ~x+y

>

> # Define 10 x 10 subplots

> grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

> survey$subplot <- over(survey,grid10)

> # B. Now find fraction of tree crown overlapping each subplot

> for (i in 1:100) {

> # Extract centroïd of each the ith cell

> centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

> corner <-

> data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

>

>

> # Find trees in a max radius (define above)

> tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]

>

>

> # Define tree crown based on tree diameter

> tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

>

> # Compute the distance from each tree to cell's borders

> pDist <- vector()

> for (k in 1:nrow(tem)) {

> pDist[k] <-

> gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))

> }

>

> # Keeps only trees whose crown is lower than the above

> distance (=overlap)

> overlap.trees <- tem[which(pDist<=tem$crownr),]

> overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown area

>

> # Creat polygons from overlapping crowns

> c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

> lonlat=F, dissolve=F)

> crown <- polygons(c1)

> Crown <-

> SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

>

> # Create a fine grid points to retrieve the fraction of

> overlapping crowns

> max.dist <- ceiling(sqrt(which.max((centro$x -

> overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

> to narrow search

>

> finegrid <-

> as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

> coordinates(finegrid) <- ~ x+y

> A <- extract(Crown,finegrid)

> Crown@data$ID <- seq(1,length(crown),1)

> B <- as.data.frame(table(A$poly.ID))

> B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

> B$overlap <- B$Freq/B$crown.area

> B$overlap[B$overlap>1] <- 1

> res[i] <- sum(B$overlap)

> }

> # C. Check the result

> res # sum of crown fraction overlapping each cell (works fine)

>

--

Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of

Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032

Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web:

gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

### Optimizing spatial query in R

Dear All,

I am desperately trying to fasten my algorithm to estimate the fraction of

tree crown that overlap a given 10x10 subplot in a forest plot. I have

combined a set of spatial functions (gDistance, extract) & objects

(SpatialGrid, SpatialPolygons) in a way that is probably not the most

efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

forest plot).

My detailed problem and a reproducible example are posted on Stackoverflow

<http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r> and

append below, if you wanna have a look. Apart of Amazon Web Server, is

anyone aware of a sever where to execute (and save results back) R codes

online?

Thanks for any help.

Best regards,

# A. Define objects

require(sp)

require(raster)

require(rgdal)

require(rgeos)

require(dismo)

radius=25 # max search radius around 10 x 10 m cells

res <- vector() # where to store results

# Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)

set.seed(0)

survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

coordinates(survey) <- ~x+y

# Define 10 x 10 subplots

grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

survey$subplot <- over(survey,grid10)

# B. Now find fraction of tree crown overlapping each subplot

for (i in 1:100) {

# Extract centroïd of each the ith cell

centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

corner <-

data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

# Find trees in a max radius (define above)

tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]

# Define tree crown based on tree diameter

tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

# Compute the distance from each tree to cell's borders

pDist <- vector()

for (k in 1:nrow(tem)) {

pDist[k] <-

gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))

}

# Keeps only trees whose crown is lower than the above

distance (=overlap)

overlap.trees <- tem[which(pDist<=tem$crownr),]

overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown area

# Creat polygons from overlapping crowns

c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

lonlat=F, dissolve=F)

crown <- polygons(c1)

Crown <-

SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

# Create a fine grid points to retrieve the fraction of

overlapping crowns

max.dist <- ceiling(sqrt(which.max((centro$x -

overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

to narrow search

finegrid <-

as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

coordinates(finegrid) <- ~ x+y

A <- extract(Crown,finegrid)

Crown@data$ID <- seq(1,length(crown),1)

B <- as.data.frame(table(A$poly.ID))

B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

B$overlap <- B$Freq/B$crown.area

B$overlap[B$overlap>1] <- 1

res[i] <- sum(B$overlap)

}

# C. Check the result

res # sum of crown fraction overlapping each cell (works fine)

--

_____________________________

Ervan Rutishauser, PhD

<http://carboforexpert.ch/>STRI post-doctoral fellow

CarboForExpert

<http://carboforexpert.ch/>

Skype: ervan-CH

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I am desperately trying to fasten my algorithm to estimate the fraction of

tree crown that overlap a given 10x10 subplot in a forest plot. I have

combined a set of spatial functions (gDistance, extract) & objects

(SpatialGrid, SpatialPolygons) in a way that is probably not the most

efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha

forest plot).

My detailed problem and a reproducible example are posted on Stackoverflow

<http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r> and

append below, if you wanna have a look. Apart of Amazon Web Server, is

anyone aware of a sever where to execute (and save results back) R codes

online?

Thanks for any help.

Best regards,

# A. Define objects

require(sp)

require(raster)

require(rgdal)

require(rgeos)

require(dismo)

radius=25 # max search radius around 10 x 10 m cells

res <- vector() # where to store results

# Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)

set.seed(0)

survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))

coordinates(survey) <- ~x+y

# Define 10 x 10 subplots

grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))

survey$subplot <- over(survey,grid10)

# B. Now find fraction of tree crown overlapping each subplot

for (i in 1:100) {

# Extract centroïd of each the ith cell

centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]

corner <-

data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

# Find trees in a max radius (define above)

tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]

# Define tree crown based on tree diameter

tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

# Compute the distance from each tree to cell's borders

pDist <- vector()

for (k in 1:nrow(tem)) {

pDist[k] <-

gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))

}

# Keeps only trees whose crown is lower than the above

distance (=overlap)

overlap.trees <- tem[which(pDist<=tem$crownr),]

overlap.trees$crowna <-overlap.trees$crownr^2*pi # compute crown area

# Creat polygons from overlapping crowns

c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,

lonlat=F, dissolve=F)

crown <- polygons(c1)

Crown <-

SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

# Create a fine grid points to retrieve the fraction of

overlapping crowns

max.dist <- ceiling(sqrt(which.max((centro$x -

overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance

to narrow search

finegrid <-

as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))

coordinates(finegrid) <- ~ x+y

A <- extract(Crown,finegrid)

Crown@data$ID <- seq(1,length(crown),1)

B <- as.data.frame(table(A$poly.ID))

B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)

B$overlap <- B$Freq/B$crown.area

B$overlap[B$overlap>1] <- 1

res[i] <- sum(B$overlap)

}

# C. Check the result

res # sum of crown fraction overlapping each cell (works fine)

--

_____________________________

Ervan Rutishauser, PhD

<http://carboforexpert.ch/>STRI post-doctoral fellow

CarboForExpert

<http://carboforexpert.ch/>

Skype: ervan-CH

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo