Subscribe to R-sig-geo feed
This is an archive for 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: 54 min 4 sec ago

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

Sun, 02/26/2017 - 04:41
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

raster brick from netcdf - specify range

Sat, 02/25/2017 - 13:03
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

Re: split SLDF with attributes

Fri, 02/24/2017 - 15:13
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

Re: split SLDF with attributes

Fri, 02/24/2017 - 09:31
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

Re: split SLDF with attributes

Fri, 02/24/2017 - 00:45
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

Re: calculate selected area from stacked raster object

Thu, 02/23/2017 - 16:44
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

split SLDF with attributes

Thu, 02/23/2017 - 13:27
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

Origin and Destination-specific doubly constrained model

Wed, 02/22/2017 - 11:03
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

calculate selected area from stacked raster object

Wed, 02/22/2017 - 10:45
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

Re: Subset a SpLineDF

Wed, 02/22/2017 - 08:12
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
signature.asc (484 bytes) Download Attachment

Subset a SpLineDF

Wed, 02/22/2017 - 04:49
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

Re: Optimizing spatial query in R

Mon, 02/20/2017 - 21:15
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

Re: Optimizing spatial query in R

Mon, 02/20/2017 - 04:04
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

Re: Optimizing spatial query in R

Mon, 02/20/2017 - 04:03
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

Re: Optimizing spatial query in R

Mon, 02/20/2017 - 03:56
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

Optimizing spatial query in R

Mon, 02/20/2017 - 02:45
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

Re: multiple lines

Sat, 02/18/2017 - 19:43
There's a number of questions I have, and a few notes.

You want line segments (two points, start,end) in x,y,time, identified by
boat?

Spatial is not really suited to this kind of multidimensional line data,
but there are specialized tools like spacetime and spatstat, no real
coherence between them but you might find a sub-domain that fits for you.
Explore the Spatiotemporal CRAN Task View for options.

In Spatial, or in sf, this would be a data set of line segment features,
with time start-end and other grouping observations, like boat ID.

There's no need to project your raster, which is costly and lossy, you can
transform the boat data to the CRS of the raster.

I may be able to help more if you clarify a bit.

Cheers, Mike

On Sat, Feb 18, 2017, 03:21 marta azores <[hidden email]> wrote:

> Hello
> I am trying to create multiple lines based on boat positions and calculate
> the distance for each boat track. I have already managed to create a single
> line with all the positions and calculate the distance but I can't managed
> to divided the line by transects which are previously defined based on the
> boat positions. These transects are not consecutive as they can be from
> different days and boats. This is an example of my data:
>
>  dput(effort)
>  structure(list(Ref = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 15L,
>  16L), Effort = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L),
>   Trip = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
>   2L), .Label = c("P001.01", "P0010.01"), class = "factor"),
>  Code = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
>  1L), .Label = "P1968", class = "factor"), Date = structure(c(1L,
>  1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("04/05/2010",
>  "19/05/2010"), class = "factor"), Time = structure(c(5L,
>  6L, 7L, 9L, 10L, 1L, 2L, 3L, 4L, 8L, 11L), .Label = c("08:12",
>  "08:25", "08:43", "08:59", "17:30", "17:51", "17:57", "18:02",
>  "18:09", "18:31", "18:39"), class = "factor"), LAT = c(38.67173,
>  38.67035, 38.66965, 38.66908, 38.64717, 38.22252, 38.26693,
>  38.18948, 38.17187, 38.25057, 38.18335), LONG = c(-28.7492,
> -28.64475, -28.70907, -28.72917, -28.63157, -28.90358, -28.87397,
> -29.07098, -28.9542, -28.8876, -28.92818), transect = c(1L, 1L, 1L,
> 1L, 1L, 2L, 2L, 2L, 2L, 4L, 4L)), .Names = c("Ref", "Effort", "Trip",
> "Code", "Date", "Time", "LAT", "LONG", "transect"), class =
> "data.frame", row.names = c(NA, -11L))
>
> The code that I have applied so far is this:
>
> ### boat points the locations of the boats are sort by date and time!
> effort <- read.table(paste0(path_data,"effort_test.csv"), header=TRUE,
>    sep=",", na.strings="NA", dec=".", strip.white=TRUE)
> pos<-as.data.frame(cbind(effort$LAT,effort$LONG, effort$transect))
> str(pos)
> sp::coordinates(pos) <- ~V2+V1
> sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
> ### project it to somewhat equal area
> effortSP <- sp::spTransform(pos, CRS("+proj=laea +lat_0=30 +lon_0=-30
> +x_0=0   +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))# plot boat
> positions
> # cost raster sea and island
> azoTS1<- raster::raster(paste0(path_data,"azoTS1.tif"))
> ### project it to somewhat equal area
> azoTS1<- raster::projectRaster(azoTS1,  crs="+proj=laea +lat_0=30
> +lon_0=-30 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
> ###scale  real cost file
> azoTS1[azoTS1 ==2]<-100
> plot(azoTS1)
> ### increase the spatial resolution
> costraw<-raster::disaggregate(azoTS1,fact=c(5, 5))
> # prepare it for analysis
> transition<- gdistance::transition(1/costraw, mean, directions=8)
> trCostS16 <- gdistance::geoCorrection(transition, type="c")
> # calculating the first segment of the whole sailing path
> line<- gdistance::shortestPath(trCostS16,
> effortSP@coords[1,],effortSP@coords[2,], output="SpatialLines")
> lines(line)
> ### here we start with the for-loopfor (i in (seq(2,length(effortSP) -
> 1))) {# calculation of the rest of the segements
> nextSegment<- gdistance::shortestPath(trCostS16,
> effortSP@coords[i,],effortSP@coords[i+1,], output="SpatialLines")#
> simple addition combines the single spatialline segements
> line <- nextSegment + line# we plot each new segement
> lines(nextSegment)}
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Calculating distances on terrain with topographic variables

Fri, 02/17/2017 - 18:41
Hello,

I'm trying to calculate pairwise distances among points, having in hand
five topographic variables: mean elevation, slope, convexity, aspect and
stage. It's a 280x220m plot divided into 20x20m quadrats (11 rows and 14
columns). So the total amount of quadrats is 154. I have these topographic
variables from each quadrat or more precisely from each x,y point at the
centre of each quadrat (x and y are simply the distance in meters from the
origin). Since it's a study of animal movement in a very steep terrain, I'd
like to get the actual distance between each pair of quadrats (x,y points)
and an index of terrain resistance or something like that.

I've tried different modeling methods, using gdistance and its least-cost
distance analysis, for instance, but I didn't get any closer to a
satisfactory result.

I really appreciate any help,


Thanks a lot,
Dim

        [[alternative HTML version deleted]]

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

Pages