Re: [R-sig-Geo] "g.list" in rgrass7 returns an integer and not a char. vector

[Mathieu Rajerison]

Thanks a lot !

Le ven. 12 avr. 2019 à 15:05, Micha Silver  a écrit :

>
> On 12/04/2019 15:15, Mathieu Rajerison wrote:
>
> Hi List,
>
> I try to get the list of rasters which are in my working environment but
> the result of execGRASS("g.list") is an integer, not a character vector
>
> Why is that ?
>
>
> execGRASS("g.list", type="raster", pattern="global*") %>% class()
>
> global_glob_rad_day1
> global_glob_rad_day2
> global_glob_rad_day3
> [1] "integer"
>
> How to get a list of my rasters as a character ?
>
>
> Adding intern=TRUE to the execGRASS command works for me:
>
>
> > execGRASS("g.list", type="rast", pattern="meron*", intern=TRUE) %>% class
> [1] "character"
>
>
> --
> Micha Silver
> Ben Gurion Univ.
> Sde Boker, Remote Sensing Lab
> cell: +972-523-665918
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] "g.list" in rgrass7 returns an integer and not a char. vector

[Mathieu Rajerison]

Hi List,

I try to get the list of rasters which are in my working environment but
the result of execGRASS("g.list") is an integer, not a character vector

Why is that ?

> execGRASS("g.list", type="raster", pattern="global*") %>% class()
global_glob_rad_day1
global_glob_rad_day2
global_glob_rad_day3
[1] "integer"

How to get a list of my rasters as a character ?

Thaks in advance for any guidance,

Mathieu

here is my sessionInfo() :
> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
 [1] bindrcpp_0.2.2 rgrass7_0.1-12 XML_3.98-1.16
gdalUtils_2.0.1.14 forcats_0.3.0  stringr_1.3.1  dplyr_0.7.6
 [8] purrr_0.2.5readr_1.1.1tidyr_0.8.1tibble_1.4.2
 ggplot2_3.0.0  tidyverse_1.2.1raster_2.6-7
[15] sp_1.3-1   sf_0.6-3

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17  lubridate_1.7.4   lattice_0.20-35   class_7.3-14
  assertthat_0.2.0  psych_1.8.4   foreach_1.4.4
 [8] R6_2.2.2  cellranger_1.1.0  plyr_1.8.4e1071_1.6-8
 httr_1.3.1pillar_1.2.3  rlang_0.2.1
[15] lazyeval_0.2.1readxl_1.1.0  rstudioapi_0.7R.utils_2.6.0
 R.oo_1.22.0   rgdal_1.3-3   foreign_0.8-70
[22] munsell_0.5.0 broom_0.4.5   compiler_3.5.1modelr_0.1.2
  pkgconfig_2.0.1   mnormt_1.5-5  tidyselect_0.2.4
[29] codetools_0.2-15  crayon_1.3.4  withr_2.1.2
 R.methodsS3_1.7.1 grid_3.5.1nlme_3.1-137  spData_0.2.9.0
[36] jsonlite_1.5  gtable_0.2.0  DBI_1.0.0 magrittr_1.5
  units_0.6-0   scales_0.5.0  cli_1.0.0
[43] stringi_1.1.7 reshape2_1.4.3xml2_1.2.0
iterators_1.0.10  tools_3.5.1   glue_1.2.0hms_0.4.2
[50] parallel_3.5.1yaml_2.1.19   colorspace_1.3-2  classInt_0.2-3
  rvest_0.3.2   bindr_0.1.1   haven_1.1.2

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spplot : different legends per plot & mode of classification by default ?

[Mathieu Rajerison]

Thanks a lot, it's perfect
ᐧ

2016-10-03 13:22 GMT+02:00 A. Marcia BARBOSA :

> Hi,
>
> Chek out the code under fig12.R at http://rspatial.r-forge.r-proj
> ect.org/gallery/.
>
> If you don't necessarily have to use spplot, I find the choroLayer
> function in the cartography package much simpler. It uses the plot
> function, so it's compatible with par (allowing you to set mfrow, mar,
> etc.) and with other related functions (e.g. add, title), and you can
> define the 'breaks' for each map.
>
> Cheers,
> AMB
>
>
> 2016-10-03 11:18 GMT+01:00 Mathieu Rajerison 
> :
>
>> Hi,
>>
>>
>> I have several variables of numeric data that I whish to spplot.
>>
>> How to push one classification per plot ?
>>
>> Also, what is the mode of classification used by default by spplot ? Is
>> the
>> data split into quantiles ?
>>
>>
>> Thanks in advance for your answer
>>
>> Mathieu
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] DEM plot3D and overlay an aerial image

[Mathieu Rajerison]

Thanks Barry.

The doc says that drape uses a RasterLayer but not a RasterStack.

So I don't know to overlay my RGB 3-dim layer on top of the DEM.

Maybe I should convert the RGB to a unique value ?
How to convert a RGB value to a unique value and does plot3D accept it ?
ᐧ

2016-10-03 13:14 GMT+02:00 Barry Rowlingson :

> On Mon, Oct 3, 2016 at 10:38 AM, Mathieu Rajerison
>  wrote:
> > Hi R-List,
> >
> >
> > I have a DEM on one hand, and on the other hand, I have an RGB aerial
> image
> >
> > I tried rasterVis and plot3D function, but I didn't find how to use the
> > colors of my RGB aerial image.
>
> from the help for plot3D:
>
>   drape: RasterLayer, to 'drape' colors representing the values of
>   this layer on the 3D representation of layer ‘x’. In this
>   case‘x’ typically has elevation data
>
> So plot(x=mydem, drape=myimage) will make a 3d plot with the heights
> from `mydem` and the colours from `myimage` draped over it - if
> `myimage` is a 3-layer RGB raster then I suspect you'll get what you
> want...
>
> Barry
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] spplot : different legends per plot & mode of classification by default ?

[Mathieu Rajerison]

Hi,


I have several variables of numeric data that I whish to spplot.

How to push one classification per plot ?

Also, what is the mode of classification used by default by spplot ? Is the
data split into quantiles ?


Thanks in advance for your answer

Mathieu
ᐧ

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] DEM plot3D and overlay an aerial image

[Mathieu Rajerison]

Hi R-List,


I have a DEM on one hand, and on the other hand, I have an RGB aerial image

I tried rasterVis and plot3D function, but I didn't find how to use the
colors of my RGB aerial image.

For the moment, I used rgl instead although it doesn't seem very
appropriate for georeferenced data..

Here's the RGL code :

gruissan =
stack("../MAISON/DATA/geo/raster/bdortho_marseille/marseille.tif")

r = values(gruissan[[1]])/255
g = values(gruissan[[2]])/255
b = values(gruissan[[3]])/255
col = rgb(r, g, b)

library(rgl)
persp3d(1:nrow(mnt), 1:ncol(mnt), values(mnt), col=col)

Best,

Mathieu


ᐧ

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Minimum bounding circle from cluster of points

[Mathieu Rajerison]

There was the answer in R here :
http://www.psychologie.uni-kiel.de/de#draw-the-minimum-bounding-box

But the page doesn't seem to exist anymore..
ᐧ

2016-07-11 17:26 GMT+02:00 Tina Cormier :

> Wow! I'm humbled that you all took the time to help me out! You had great
> suggestions and sent along some very handy functions. I have the idea now
> of what I need to do - find the longest distance between two points in the
> cluster, the midpoint of that line is the center of the circle, and there
> are numerous ways to create the circle, but once you have the radius, you
> can just buffer the center point. Fantastically simple. Someone also
> suggested using PostGIS (ST_MinimumBoundingCircle), which I will definitely
> try. I have a postgreSQL db that I use for some simple stuff, and I'm not a
> guru by any means, but this will help me gain a little more experience.
>
> Thank you again, and I can't wait to pay if forward the next time someone
> asks about a topic I know well!
>
> Cheers,
> Tina
> ***
> In SUMMARY, functions that solve this problem are listed below:
>
> 1. Using a *3-liner* from Robert Hijmans, using optimization:
> library(raster)
> library(rgeos)
>
> set.seed(7)
> n <- 4
> xy <- cbind(runif(n), runif(n))
>
>
>
>
> *f <- function(p) { max(pointDistance(rbind(p), xy, lonlat=FALSE)) }p <-
> optim(colMeans(xy), f)cc <- buffer(SpatialPoints(rbind(p$par)),
> width=p$value, quadsegs=45)*
> plot(cc)
> points(xy, pch=20, cex=1.5)
> points(rbind(p$par), pch='*', col='red', cex=3)
> lines(rbind(p$par, c(p$par[1]+p$value, p$par[2])), col='blue', lwd=2,
> lty=2)
> text(p$par[1]+0.5*p$value, p$par[2], paste("radius =", round(p$value,
> 2)), pos=3, cex=.85)
> text(p$par[1], p$par[2], labels=paste0("(",paste(round(p$par, 2),
> collapse=", "), ")"), pos=1, cex=.75)
> ***
>
> 2. A function from Adrian Baddeley using spatstat
> library(spatstat)
>
> circumcircle <- function(x, ...) {   UseMethod("circumcircle") }
>
> circumcircle.owin <- function(x, ...) {
>   d2 <- fardist(x, ..., squared=TRUE)
>   z <- where.min(d2)
>   r <- sqrt(min(d2))
>   w <- disc(centre=z, radius=r)
>   return(w)
> }
>
> circumcircle.ppp <- function(x, ...) {
>   circumcircle(convexhull(x))
> }
> ***
>
> 3. A handy link that walks through the process in a bit longer way around,
> but still works:
> http://www.uni-kiel.de/psychologie/rexrepos/posts/diagBounding.html
>
>
> On Fri, Jul 8, 2016 at 10:31 AM, Tina Cormier 
> wrote:
>
> > Hi all,
> >
> > Looking for help aggregating some field data subplots to the plot level
> > (data currently in shapefile format). I have clusters of 4 points (4
> > subplots = 1 plot). I'd like to create a circle around each cluster that
> is
> > the smallest circle that would encompass all 4 points. Sort of like a
> > convex hull, but a circle (convex hull, in this case, would give me a
> > triangle because of the layout of the subplots). I am familiar with (and
> > use frequently) the typical geo packages in R - with rgdal, sp, raster,
> and
> > maptools being my most frequent flyers. Is there a function I'm missing
> in
> > these (or some other) packages that might be able to help me out? In the
> > attribute table, I have a unique ID for each plot/cluster. So for each
> plot
> > ID (in this case, consisting of 4 subplots), I'd like to build a circle
> > around all of the subplots.
> >
> > > str(pts)
> > Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
> >   ..@ data   :'data.frame': 160 obs. of  12 variables:
> >   .. ..$ Plot  : num [1:160] 1 1 1 1 2 2 2 2 3 3 ...
> >   .. ..$ Subplot   : num [1:160] 1 2 3 4 1 2 3 4 1 2 ...
> >
> > I should also mention that it's not always 4 points (subplots), and they
> > aren't always covering the same size area on the ground, so buffering by
> a
> > constant distance isn't the answer. I have subset out 3 plots (so 12
> > subplots) into a test file here:
> > https://dl.dropboxusercontent.com/u/72421241/test_subplotsToPlots.zip
> >
> > I have R version 3.3.0. I know it's customary to include code for what
> > I've tried, bt, that would be a blank canvas at this point. Googling
> > hasn't really been fruitful for this issue, so I thought you folks might
> > have some good ideas!
> >
> > Thanks,
> > Tina
> >
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] igraph and spatial

[Mathieu Rajerison]

In this example, I generate a bezier line with random curvature within a
specific range :

from = mst.bh[i, 1] # coordinates of the source node
to = mst.bh[i, 2] # coordinates of the target node

d = sqrt((coords[from, 1] - coords[to, 1])^2 + (coords[from, 2] -
coords[to, 2])^2) # the distance bewteen the two nodes

midpoint = c((coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) >
0, 1, -1) * runif(1,0,d/3),
 (coords[from, 2] + coords[to, 2])/2 + ifelse(runif(1,-1,1)
> 0, 1, -1) * runif(1,0,d/3))

 b = bezier(t=t <- seq(0, 1, length=100), p=rbind(coords[from, ], midpoint,
coords[to, ]))

rowname = paste(from, to, sep="_")
l = makeLineFromCoords(b, rowname)

2016-02-25 11:55 GMT+01:00 Mathieu Rajerison :

> For curved lines, you could use the bezier library.
>
> In this example, I generate a bezier line with random curvature within a
> specific range :
>
> from = mst.bh[i, 1] # coordinates of the source node
> to = mst.bh[i, 2] # coordinates of the target node
>
> d = sqrt((coords[from, 1] - coords[to, 1])^2 + (coords[from, 2] -
> coords[to, 2])^2) # the distance bewteen the two nodes
>
> midpoint = c((coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) >
> 0, 1, -1) * runif(1,0,d/3),
>  (coords[from, 2] + coords[to, 2])/2 +
> ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3))
>
>  b = bezier(t=t <- seq(0, 1, length=100), p=rbind(coords[from, ],
> midpoint, coords[to, ]))
>
> rowname = paste(from, to, sep="_")
> l = makeLineFromCoords(b, rowname)
>
> 2016-02-25 11:44 GMT+01:00 Luca Candeloro :
>
>> Thank you Mathieu,
>> in the way you propose, however, you get straight lines...
>> I'd like instead to have curved edges.
>> any suggestions?
>> cheers, luca
>>
>> 2016-02-25 11:39 GMT+01:00 Mathieu Rajerison > >:
>>
>>> Hi,
>>>
>>> Recently, I used igraph with spatial objects.
>>>
>>> Once you have the coordinates of the nodes composing an edge, you can
>>> transform it into a spatiallines object with this function :
>>>
>>> makeLineFromCoords <- function(coords, i) {
>>>   Sl1 = Line(coords)
>>>   S1 = Lines(list(Sl1), ID=as.character(i))
>>>   Sl = SpatialLines(list(S1))
>>>   return(Sl)
>>> }
>>>
>>> For all the routing stuff within R, this resource by B. Rowlingson is
>>> quite inspiring :
>>>
>>> http://rstudio-pubs-static.s3.amazonaws.com/1572_7599552b60454033a0d5c5e6d2e34ffb.html
>>>
>>> best,
>>>
>>> Mathieu
>>>
>>> 2016-02-22 9:30 GMT+01:00 Michael Sumner :
>>>
>>>> On Mon, 22 Feb 2016 at 19:13 Luca Candeloro 
>>>> wrote:
>>>>
>>>> > Dear all,
>>>> > plotting graphs with given vertex coordinates is possible inside
>>>> igraph
>>>> > plot function.
>>>> > I would like to transform vilsualized edges into
>>>> SpatialLinesDataFrame, so
>>>> > that they could be exported in .shp file...
>>>> > Any suggestion?
>>>> >
>>>>
>>>> Can you please provide a working igraph example that we can start with?
>>>>
>>>> Do you want coordinates as seen in a particular plot, or can it be more
>>>> arbitrary?
>>>>
>>>> Cheers, Mike.
>>>>
>>>>
>>>>
>>>> > thanks in advance,
>>>> > Luca
>>>> >
>>>> > [[alternative HTML version deleted]]
>>>> >
>>>> > ___
>>>> > R-sig-Geo mailing list
>>>> > R-sig-Geo@r-project.org
>>>> > 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
>>>> R-sig-Geo@r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] igraph and spatial

[Mathieu Rajerison]

For curved lines, you could use the bezier library.

In this example, I generate a bezier line with random curvature within a
specific range :

from = mst.bh[i, 1] # coordinates of the source node
to = mst.bh[i, 2] # coordinates of the target node

d = sqrt((coords[from, 1] - coords[to, 1])^2 + (coords[from, 2] -
coords[to, 2])^2) # the distance bewteen the two nodes

midpoint = c((coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) >
0, 1, -1) * runif(1,0,d/3),
 (coords[from, 2] + coords[to, 2])/2 + ifelse(runif(1,-1,1)
> 0, 1, -1) * runif(1,0,d/3))

 b = bezier(t=t <- seq(0, 1, length=100), p=rbind(coords[from, ], midpoint,
coords[to, ]))

rowname = paste(from, to, sep="_")
l = makeLineFromCoords(b, rowname)

2016-02-25 11:44 GMT+01:00 Luca Candeloro :

> Thank you Mathieu,
> in the way you propose, however, you get straight lines...
> I'd like instead to have curved edges.
> any suggestions?
> cheers, luca
>
> 2016-02-25 11:39 GMT+01:00 Mathieu Rajerison 
> :
>
>> Hi,
>>
>> Recently, I used igraph with spatial objects.
>>
>> Once you have the coordinates of the nodes composing an edge, you can
>> transform it into a spatiallines object with this function :
>>
>> makeLineFromCoords <- function(coords, i) {
>>   Sl1 = Line(coords)
>>   S1 = Lines(list(Sl1), ID=as.character(i))
>>   Sl = SpatialLines(list(S1))
>>   return(Sl)
>> }
>>
>> For all the routing stuff within R, this resource by B. Rowlingson is
>> quite inspiring :
>>
>> http://rstudio-pubs-static.s3.amazonaws.com/1572_7599552b60454033a0d5c5e6d2e34ffb.html
>>
>> best,
>>
>> Mathieu
>>
>> 2016-02-22 9:30 GMT+01:00 Michael Sumner :
>>
>>> On Mon, 22 Feb 2016 at 19:13 Luca Candeloro 
>>> wrote:
>>>
>>> > Dear all,
>>> > plotting graphs with given vertex coordinates is possible inside igraph
>>> > plot function.
>>> > I would like to transform vilsualized edges into
>>> SpatialLinesDataFrame, so
>>> > that they could be exported in .shp file...
>>> > Any suggestion?
>>> >
>>>
>>> Can you please provide a working igraph example that we can start with?
>>>
>>> Do you want coordinates as seen in a particular plot, or can it be more
>>> arbitrary?
>>>
>>> Cheers, Mike.
>>>
>>>
>>>
>>> > thanks in advance,
>>> > Luca
>>> >
>>> > [[alternative HTML version deleted]]
>>> >
>>> > ___
>>> > R-sig-Geo mailing list
>>> > R-sig-Geo@r-project.org
>>> > 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
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] igraph and spatial

[Mathieu Rajerison]

Hi,

Recently, I used igraph with spatial objects.

Once you have the coordinates of the nodes composing an edge, you can
transform it into a spatiallines object with this function :

makeLineFromCoords <- function(coords, i) {
  Sl1 = Line(coords)
  S1 = Lines(list(Sl1), ID=as.character(i))
  Sl = SpatialLines(list(S1))
  return(Sl)
}

For all the routing stuff within R, this resource by B. Rowlingson is quite
inspiring :
http://rstudio-pubs-static.s3.amazonaws.com/1572_7599552b60454033a0d5c5e6d2e34ffb.html

best,

Mathieu

2016-02-22 9:30 GMT+01:00 Michael Sumner :

> On Mon, 22 Feb 2016 at 19:13 Luca Candeloro 
> wrote:
>
> > Dear all,
> > plotting graphs with given vertex coordinates is possible inside igraph
> > plot function.
> > I would like to transform vilsualized edges into SpatialLinesDataFrame,
> so
> > that they could be exported in .shp file...
> > Any suggestion?
> >
>
> Can you please provide a working igraph example that we can start with?
>
> Do you want coordinates as seen in a particular plot, or can it be more
> arbitrary?
>
> Cheers, Mike.
>
>
>
> > thanks in advance,
> > Luca
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > 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
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] RStoolbox loading, Rscript, printing and e acute (é)

[Mathieu Rajerison]

I forgot to tell it's only when calling the script from Rscript that the e
acute isn't printed correctly.

in my case, I'm under Windows

R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C   LC_TIME=French_France.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] RStoolbox_0.1.1 raster_2.3-40   sp_1.0-17

loaded via a namespace (and not attached):
 [1] car_2.0-25 caret_6.0-57   codetools_0.2-9
 colorspace_1.2-6   digest_0.6.8   doParallel_1.0.8
 [7] foreach_1.4.2  geosphere_1.4-3ggplot2_1.0.1  grid_3.1.2
  gtable_0.1.2   iterators_1.0.7
[13] lattice_0.20-29lme4_1.1-10magrittr_1.5   MASS_7.3-35
 Matrix_1.2-2   MatrixModels_0.4-1
[19] mgcv_1.8-3 minqa_1.2.4munsell_0.4.2  nlme_3.1-118
  nloptr_1.0.4   nnet_7.3-8
[25] parallel_3.1.2 pbkrtest_0.4-2 plyr_1.8.3 proto_0.3-10
  quantreg_5.19  Rcpp_0.12.0
[31] reshape2_1.4.1 rgeos_0.3-8scales_0.2.5   SparseM_1.7
 splines_3.1.2  stats4_3.1.2
[37] stringi_0.5-5  stringr_1.0.0  tools_3.1.2XML_3.98-1.3


2016-01-21 14:55 GMT+01:00 Edzer Pebesma :

> Mathieu, I can't reproduce this on ubuntu 14.04:
>
> >
> > cat("first \U00E9 \n")
> first é
> > cat("--before--\n")
> --before--
> > cat(paste("encoding = ", getOption("encoding"), "\n"))
> encoding =  native.enc
> > cat(paste("locale = ", Sys.getlocale(), "\n"))
> locale =
>
> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_GB.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_GB.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_GB.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_GB.UTF-8;LC_IDENTIFICATION=C
>
> > # suppressMessages(suppressWarnings(require("RStoolbox", quietly =
> TRUE)))
> > # provoque erreurs accent e aigu
> > require("RStoolbox")
> Loading required package: RStoolbox
> > cat("\n")
>
> > cat("--after--\n")
> --after--
> > cat(paste("encoding = ", getOption("encoding"),"\n"))
> encoding =  native.enc
> > cat(paste("locale = ", Sys.getlocale(), "\n"))
> locale =
>
> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_GB.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_GB.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_GB.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_GB.UTF-8;LC_IDENTIFICATION=C
>
> > cat(" second \U00E9")
>  second é> sessionInfo()
> R version 3.2.3 (2015-12-10)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Ubuntu 14.04.3 LTS
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
>  [3] LC_TIME=en_GB.UTF-8LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_GB.UTF-8LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_GB.UTF-8   LC_NAME=C
>  [9] LC_ADDRESS=C   LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
>
> other attached packages:
> [1] RStoolbox_0.1.3
>
> loaded via a namespace (and not attached):
>  [1] Rcpp_0.12.2raster_2.5-2   magrittr_1.5
> splines_3.2.3
>  [5] MASS_7.3-44doParallel_1.0.10  munsell_0.4.2
> geosphere_1.5-1
>  [9] colorspace_1.2-6   lattice_0.20-33foreach_1.4.3
> minqa_1.2.4
> [13] stringr_1.0.0  car_2.1-0  plyr_1.8.3
> tools_3.2.3
> [17] parallel_3.2.3 nnet_7.3-11pbkrtest_0.4-2
> caret_6.0-64
> [21] grid_3.2.3 gtable_0.1.2   nlme_3.1-123   mgcv_1.8-7
>
> [25] quantreg_5.19  rgeos_0.3-15   MatrixModels_0.4-1
> iterators_1.0.8
> [29] lme4_1.1-10Matrix_1.2-3   nloptr_1.0.4
> reshape2_1.4.1
> [33] ggplot2_2.0.0  codetools_0.2-14   sp_1.2-2
> stringi_1.0-1
> [37] scales_0.3.0   XML_3.98-1.3   stats4_3.2.3   SparseM_1.7
>
> On 21/01/16 13:58, Mathieu Rajerison wrote:
> > Hi,
> >
> > This post is not spatially-oriented but deals with an excellent spatial R
> > package named RStoobox for remote sensing.
> >
> > I have a process that prints advancement operations with cat() functions.
> > Some contain e acute as it is for french administrations.
> >
> > I noticed I had problems with the writing of e acutes after loading
> > RStoolbox package. I noticed the same with igraph as well
> >
> > If I write an R file with the following lines :
>

[R-sig-Geo] RStoolbox loading, Rscript, printing and e acute (é)

[Mathieu Rajerison]

Hi,

This post is not spatially-oriented but deals with an excellent spatial R
package named RStoobox for remote sensing.

I have a process that prints advancement operations with cat() functions.
Some contain e acute as it is for french administrations.

I noticed I had problems with the writing of e acutes after loading
RStoolbox package. I noticed the same with igraph as well

If I write an R file with the following lines :

cat("first \U00E9")
require("RStoolbox")
cat(" second \U00E9")

It gives

first é
second Ú

If I load another package like raster, the problem doesn't appear.

To see if there was any change relative to encodings, I wrote an R file
with the following lines which gives the system encoding before and after :

cat("first \U00E9 \n")
cat("--before--\n")
cat(paste("encoding = ", getOption("encoding"), "\n"))
cat(paste("locale = ", Sys.getlocale(), "\n"))
# suppressMessages(suppressWarnings(require("RStoolbox", quietly = TRUE)))
# provoque erreurs accent e aigu
require("RStoolbox")
cat("\n")
cat("--after--\n")
cat(paste("encoding = ", getOption("encoding"),"\n"))
cat(paste("locale = ", Sys.getlocale(), "\n"))
cat(" second \U00E9")

But there isn't any change :

first é
--before--
encoding =  native.enc
locale =
 LC_COLLATE=French_France.1252;LC_CTYPE=French_France.1252;LC_MONETARY=
French_France.1252;LC_NUMERIC=C;LC_TIME=French_France.1252
Le chargement a nécessité le package : RStoolbox
Message d'avis :
le package 'RStoolbox' a ÚtÚ compilÚ avec la version R 3.1.3

--after--
encoding =  native.enc
locale =
 LC_COLLATE=French_France.1252;LC_CTYPE=French_France.1252;LC_MONETARY=
French_France.1252;LC_NUMERIC=C;LC_TIME=French_France.1252
 second Ú

So, i wonder where does this come from ? And how to deal with it

Thanks in advance for your answers,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions

[Mathieu Rajerison]

Hi,

And thanks for all the answers.

I've tested Robert's elegant solution successfully.
I haven't yet tested your solution, Banjamin, but it seems a really good
improvement for speeding all this up and I'll soon give it a try. Thanks
for your awesome RSToolbox package

To answer your question "why setting all pixels to 1", it's just that
afterwards, I smooth the result with a mean-shift focal function. This way,
it gives all the extracted pixels the same value/weight.

Thanks again and happy new yeaR

Mathieu


2015-12-25 6:46 GMT+01:00 Robert J. Hijmans :

> How about this:
>
> b <- brick(system.file("external/rlogo.grd", package="raster"))
> classified <- b / 255
>
> threshs <- c(.1, .3, .5)
> x <- classified < threshs
> y <- max(x, na.rm=TRUE)
> z <- reclassify(y, cbind(0, NA))
>
>
> Robert
>
>
> On Tue, Dec 22, 2015 at 1:57 AM, Mathieu Rajerison
>  wrote:
> > Hi,
> >
> >
> > I use RSToolBox to classify a RGB raster.
> >
> > I have a resulting RasterBrick which has as many layer as end members, in
> > my case 3 for different tones of blue.
> >
> > I reclassify each band with calc to extract the pixels which have a small
> > angle mapping value. The threshold used is different depending on the
> > endmember layer.
> >
> > I finally assembly all the bands with the cover function.
> >
> > I needed to increase the memory limit assigned to R to have it worked. I
> > suspect that my code could be optimized, but I don't know in which way.
> >
> > Here is the part of my code, that I think, could be optilmized, if you
> want
> > to have a look and give some advice :
> >
> > # RECLASSIFY
> > # classified is the classified RasterStack
> > # here I change the values of each band to 1 or NA depending on the
> > spectral angle mapping value.
> > # Is calc() slower than reclassify() for this purpose as I have only one
> > threshold value ?
> >
> > threshs = c(.1,.2,.1)
> > for (i in 1:nlayers(classified )) {
> >
> > clas = classified[[i]]
> > thresh=threshs[i]
> >
> > out[[i]] = calc(clas, function(x) {x[x >= thresh] = NA;
> >x[x < thresh] = 1;
> >return(x)})
> >   }
> >
> > # COVERING
> >   r = out[[1]]
> >   for(i in 2:length(out)) {
> > r = cover(out[[i]], r)   ## I cover by iteration
> >   }
> >
> > plot(r) # r is the final combined raster
> >
> > 
> > My sessionInfo() :
> >> sessionInfo()
> > R version 3.1.2 (2014-10-31)
> > Platform: x86_64-w64-mingw32/x64 (64-bit)
> >
> > locale:
> > [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
> >  LC_MONETARY=French_France.1252
> > [4] LC_NUMERIC=C   LC_TIME=French_France.1252
> >
> > attached base packages:
> > [1] stats graphics  grDevices utils datasets  methods   base
> >
> > other attached packages:
> >  [1] R.utils_2.1.0R.oo_1.19.0  R.methodsS3_1.7.0
> >  igraph_1.0.1 scatterplot3d_0.3-36
> >  [6] gdalUtils_2.0.1.7spdep_0.5-88 Matrix_1.2-2
> > maptools_0.8-34  spgrass6_0.8-8
> > [11] XML_3.98-1.3 rgeos_0.3-8  FNN_1.1
> >  rgdal_0.9-2  RStoolbox_0.1.1
> > [16] raster_2.3-40sp_1.0-17
> >
> > loaded via a namespace (and not attached):
> >  [1] boot_1.3-13car_2.0-25 caret_6.0-57   coda_0.18-1
> >  codetools_0.2-9colorspace_1.2-6
> >  [7] deldir_0.1-9   digest_0.6.8   doParallel_1.0.8
>  foreach_1.4.2
> >  foreign_0.8-61 geosphere_1.4-3
> > [13] ggplot2_1.0.1  grid_3.1.2 gtable_0.1.2
> > iterators_1.0.7lattice_0.20-29LearnBayes_2.15
> > [19] lme4_1.1-10magrittr_1.5   MASS_7.3-35
> >  MatrixModels_0.4-1 mgcv_1.8-3 minqa_1.2.4
> > [25] munsell_0.4.2  nlme_3.1-118   nloptr_1.0.4   nnet_7.3-8
> >   parallel_3.1.2 pbkrtest_0.4-2
> > [31] plyr_1.8.3 proto_0.3-10   quantreg_5.19  Rcpp_0.12.0
> >  reshape2_1.4.1 scales_0.2.5
> > [37] SparseM_1.7splines_3.1.2  stats4_3.1.2
>  stringi_0.5-5
> >  stringr_1.0.0  tools_3.1.2
> >
> > 
> > Best,
> >
> > Mathieu
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions

[Mathieu Rajerison]

Hi Benjamin,

Thanks for your answer

Sorry for not being precise. Actually, I want to threshold each band given
thresholds. I give the 1 value to each pixel, in each band, that is below
the threshold value.

Then, I want to cover all these bands into one, in order to combine all the
layers.

I tried different functions. Here are the different performances, if you
want to have a look. Te second function below is the fastest one. If you
think you have better, do not hesitate.

## FIRST FUNCTION

  fun1 = function() {
threshs = c(.1,.2,.1)
tamis = raster(classified); tamis[] = 0
out = stack(lapply(1:nlayers(classified), function(i)
cover(clamp(classified[[i]], upper = threshs[[i]], useValues = FALSE),
tamis)))
r = calc(out, sum)
r[r[]==0] = NA; r[which(!is.na(r[]))] = 1
  }

## SECOND

  fun2 = function() {
threshs = c(.1,.2,.1)
out = lapply(1:nlayers(classified), function(i) clamp(classified[[i]],
upper = threshs[[i]], useValues = FALSE))
r = out[[1]]
for(i in 2:length(out)) {
  r = cover(out[[i]], r)
}
r[which(!is.na(r[]))] = 1
  }

## THIRD

  fun3 = function() {
threshs = c(.1,.2,.1)
out=lapply(1:nlayers(classified), function(i) {
  out[[i]] = calc(classified[[i]], function(x) {x[x >= threshs[i]] = NA;
x[x < threshs[i]] = 1;
return(x)})
})
r = out[[1]]
for(i in 2:length(out)) {
  r = cover(out[[i]], r)
}
  }

  system.time(fun1())
  system.time(fun2())
  system.time(fun3())

> system.time(fun1())
   user  system elapsed
  63.052.67   65.81
>   system.time(fun2())
   user  system elapsed
  43.140.88  * 44.52 *
>   system.time(fun3())
   user  system elapsed
  48.671.51   50.62


Best,

Mathieu



2015-12-22 14:09 GMT+01:00 Benjamin Leutner <
benjamin.leut...@uni-wuerzburg.de>:

> Hi Mathieu,
>
> your question is rather difficult to understand. From the context I gather
> that you are referring to the results of the sam() function from RStoolbox.
> Further, I assume you want to threshold each layer for a maximum spectral
> angle  and then find the class with the minimum spectral angle per pixel,
> right?
>
> In this case you could do:
>
> out   <- stack(lapply(1:nlayers(classified), function(i)
> clamp(classified[[i]], upper = threshs[[i]], useValues = FALSE)))
> class <- which.min(out)
>
> Cheers,
> Benjamin
>
>
> On 22.12.2015 10:57, Mathieu Rajerison wrote:
>
>> Hi,
>>
>>
>> I use RSToolBox to classify a RGB raster.
>>
>> I have a resulting RasterBrick which has as many layer as end members, in
>> my case 3 for different tones of blue.
>>
>> I reclassify each band with calc to extract the pixels which have a small
>> angle mapping value. The threshold used is different depending on the
>> endmember layer.
>>
>> I finally assembly all the bands with the cover function.
>>
>> I needed to increase the memory limit assigned to R to have it worked. I
>> suspect that my code could be optimized, but I don't know in which way.
>>
>> Here is the part of my code, that I think, could be optilmized, if you
>> want
>> to have a look and give some advice :
>>
>> # RECLASSIFY
>> # classified is the classified RasterStack
>> # here I change the values of each band to 1 or NA depending on the
>> spectral angle mapping value.
>> # Is calc() slower than reclassify() for this purpose as I have only one
>> threshold value ?
>>
>> threshs = c(.1,.2,.1)
>> for (i in 1:nlayers(classified )) {
>>
>>  clas = classified[[i]]
>>  thresh=threshs[i]
>>
>>  out[[i]] = calc(clas, function(x) {x[x >= thresh] = NA;
>> x[x < thresh] = 1;
>> return(x)})
>>}
>>
>> # COVERING
>>r = out[[1]]
>>for(i in 2:length(out)) {
>>  r = cover(out[[i]], r)   ## I cover by iteration
>>}
>>
>> plot(r) # r is the final combined raster
>>
>> 
>> My sessionInfo() :
>>
>>> sessionInfo()
>>>
>> R version 3.1.2 (2014-10-31)
>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>>
>> locale:
>> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
>>   LC_MONETARY=French_France.1252
>> [4] LC_NUMERIC=C   LC_TIME=French_France.1252
>>
>> attached base packages:
>> [1] stats graphics  grDevices utils datasets  methods   base
>>
>> other attached packages:
>>   [1] R.utils_2.1.0R.o

[R-sig-Geo] Needing to speed up a process involving calc() and cover() raster functions

[Mathieu Rajerison]

Hi,


I use RSToolBox to classify a RGB raster.

I have a resulting RasterBrick which has as many layer as end members, in
my case 3 for different tones of blue.

I reclassify each band with calc to extract the pixels which have a small
angle mapping value. The threshold used is different depending on the
endmember layer.

I finally assembly all the bands with the cover function.

I needed to increase the memory limit assigned to R to have it worked. I
suspect that my code could be optimized, but I don't know in which way.

Here is the part of my code, that I think, could be optilmized, if you want
to have a look and give some advice :

# RECLASSIFY
# classified is the classified RasterStack
# here I change the values of each band to 1 or NA depending on the
spectral angle mapping value.
# Is calc() slower than reclassify() for this purpose as I have only one
threshold value ?

threshs = c(.1,.2,.1)
for (i in 1:nlayers(classified )) {

clas = classified[[i]]
thresh=threshs[i]

out[[i]] = calc(clas, function(x) {x[x >= thresh] = NA;
   x[x < thresh] = 1;
   return(x)})
  }

# COVERING
  r = out[[1]]
  for(i in 2:length(out)) {
r = cover(out[[i]], r)   ## I cover by iteration
  }

plot(r) # r is the final combined raster


My sessionInfo() :
> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C   LC_TIME=French_France.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
 [1] R.utils_2.1.0R.oo_1.19.0  R.methodsS3_1.7.0
 igraph_1.0.1 scatterplot3d_0.3-36
 [6] gdalUtils_2.0.1.7spdep_0.5-88 Matrix_1.2-2
maptools_0.8-34  spgrass6_0.8-8
[11] XML_3.98-1.3 rgeos_0.3-8  FNN_1.1
 rgdal_0.9-2  RStoolbox_0.1.1
[16] raster_2.3-40sp_1.0-17

loaded via a namespace (and not attached):
 [1] boot_1.3-13car_2.0-25 caret_6.0-57   coda_0.18-1
 codetools_0.2-9colorspace_1.2-6
 [7] deldir_0.1-9   digest_0.6.8   doParallel_1.0.8   foreach_1.4.2
 foreign_0.8-61 geosphere_1.4-3
[13] ggplot2_1.0.1  grid_3.1.2 gtable_0.1.2
iterators_1.0.7lattice_0.20-29LearnBayes_2.15
[19] lme4_1.1-10magrittr_1.5   MASS_7.3-35
 MatrixModels_0.4-1 mgcv_1.8-3 minqa_1.2.4
[25] munsell_0.4.2  nlme_3.1-118   nloptr_1.0.4   nnet_7.3-8
  parallel_3.1.2 pbkrtest_0.4-2
[31] plyr_1.8.3 proto_0.3-10   quantreg_5.19  Rcpp_0.12.0
 reshape2_1.4.1 scales_0.2.5
[37] SparseM_1.7splines_3.1.2  stats4_3.1.2   stringi_0.5-5
 stringr_1.0.0  tools_3.1.2


Best,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spgrass6, readVECT and overwriting ?

[Mathieu Rajerison]

Yes, you're right. I should export directly from GRASS.

When doing :

for (i in 1:3) {
  f <- readVECT("f")
}

I get :
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer  already exists in OGR data source
'D:/GRASSDB/paca/mapset/.tmp'
OGR data source with driver: ESRI Shapefile
Source: "D:/GRASSDB/paca/mapset/.tmp", layer: "f"
with 1 features
It has 5 fields
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer  already exists in OGR data source
'D:/GRASSDB/paca/mapset/.tmp'
OGR data source with driver: ESRI Shapefile
Source: "D:/GRASSDB/paca/mapset/.tmp", layer: "f"
with 1 features
It has 5 fields
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer  already exists in OGR data source
'D:/GRASSDB/paca/mapset/.tmp'
OGR data source with driver: ESRI Shapefile
Source: "D:/GRASSDB/paca/mapset/.tmp", layer: "f"
with 1 features
It has 5 fields
Warning messages:
1: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1
2: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1
3: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1

-
sessionInfo() :
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C   LC_TIME=French_France.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
 [1] R.utils_2.1.0 R.oo_1.19.0   R.methodsS3_1.7.0 igraph_1.0.1
 gdalUtils_2.0.1.7
 [6] spdep_0.5-88  Matrix_1.2-2  maptools_0.8-34   spgrass6_0.8-8
 XML_3.98-1.3
[11] rgeos_0.3-8   FNN_1.1   rgdal_0.9-2   raster_2.3-40
  sp_1.0-17
[16] RStoolbox_0.1.1

loaded via a namespace (and not attached):
 [1] boot_1.3-13car_2.0-25 caret_6.0-57   coda_0.18-1
 codetools_0.2-9
 [6] colorspace_1.2-6   deldir_0.1-9   digest_0.6.8
doParallel_1.0.8   foreach_1.4.2
[11] foreign_0.8-61 geosphere_1.4-3ggplot2_1.0.1  grid_3.1.2
  gtable_0.1.2
[16] iterators_1.0.7lattice_0.20-29LearnBayes_2.15lme4_1.1-10
 magrittr_1.5
[21] MASS_7.3-35MatrixModels_0.4-1 mgcv_1.8-3 minqa_1.2.4
 munsell_0.4.2
[26] nlme_3.1-118   nloptr_1.0.4   nnet_7.3-8
parallel_3.1.2 pbkrtest_0.4-2
[31] plyr_1.8.3 proto_0.3-10   quantreg_5.19  Rcpp_0.12.0
 reshape2_1.4.1
[36] scales_0.2.5   SparseM_1.7    splines_3.1.2  stats4_3.1.2
  stringi_0.5-5
[41] stringr_1.0.0  tools_3.1.2

2015-10-28 16:16 GMT+01:00 Rainer M Krug :

> Mathieu Rajerison  writes:
>
> > Before was the code that works.
> >
> > Below is the code that generates issues :
> >
> > (sorry the multiple posts)
> >
> > for (id in 1:length(grids)) {
> >
> > (...)
> >
> > ## SMOOTH
> >   print(">> SMOOTH")
> >   smoothed = smooth(sieved, 9)
> >   values(smoothed)[values(smoothed) > 0] = 1
> >   values(smoothed)[values(smoothed)==0] = NA
> >
> >   ## THIN
> >   # WRITE TO GRASS
> >   writeRAST(as(smoothed, "SpatialGridDataFrame"), "bdtopo", zcol="layer",
> > flags=c("overwrite"))
> >   execGRASS("g.region", rast="bdtopo")
> >   execGRASS("r.thin", parameters=list(input="bdtopo", output="thin"),
> > flags=c("overwrite"))
> >   thin = raster(readRAST("thin"))
> >
> >   ## CLEAN
> >   execGRASS("r.to.vect", parameters=list(input="thin", output="thin",
> > feature="line"), flags=c("overwrite"))
> >   execGRASS("v.clean", parameters=list(input="thin", tool="rmdangle",
> > output="clean", thresh=20), flags=c("overwrite"))
> >
> >   cleaned = readVECT("clean")
> >
> >   ## EXPORT
> >   writeOGR(cleaned, "OUT/CLEANED", id, "ESRI Shapefile")
> > }
>
> OK. So what you are doing is in a loop
>
>   write a raster to gr

Re: [R-sig-Geo] spgrass6, readVECT and overwriting ?

[Mathieu Rajerison]

Before was the code that works.

Below is the code that generates issues :

(sorry the multiple posts)

for (id in 1:length(grids)) {

(...)

## SMOOTH
  print(">> SMOOTH")
  smoothed = smooth(sieved, 9)
  values(smoothed)[values(smoothed) > 0] = 1
  values(smoothed)[values(smoothed)==0] = NA

  ## THIN
  # WRITE TO GRASS
  writeRAST(as(smoothed, "SpatialGridDataFrame"), "bdtopo", zcol="layer",
flags=c("overwrite"))
  execGRASS("g.region", rast="bdtopo")
  execGRASS("r.thin", parameters=list(input="bdtopo", output="thin"),
flags=c("overwrite"))
  thin = raster(readRAST("thin"))

  ## CLEAN
  execGRASS("r.to.vect", parameters=list(input="thin", output="thin",
feature="line"), flags=c("overwrite"))
  execGRASS("v.clean", parameters=list(input="thin", tool="rmdangle",
output="clean", thresh=20), flags=c("overwrite"))

  cleaned = readVECT("clean")

  ## EXPORT
  writeOGR(cleaned, "OUT/CLEANED", id, "ESRI Shapefile")
}

2015-10-28 14:18 GMT+01:00 Mathieu Rajerison :

> Mr Krug,
>
>
> It is difficult to reproduce an example from the spearfish dataset as it
> doesn't natively contain data in it close to the one I used.
>
> To give you th context, my goal is to detect all blue features : rivers
> from a topographic map, ie 3-band RGB raster.
>
> With RStoolbox, I perform a Spectral Angle Classification of my data using
> sample points located on (near)blue pixels
>
> With GRASS, I then thin, vectorize and clean the result.
>
> As the process must be performed on a very large scale, I split it using
> grids where sample RGB parts of the topo map are extracted with
> gdal_translate. Everything is executed within a loop. At the end of the
> loop, I export the resulted cleant shapefile in a folder. At the end, I
> have as many shapefiles as processing grids.
>
> That's why I use temporary vector datasets inside GRASS that are
> overwritten at each loop instead of storing each individual result. I
> noticed that in my process, the vector data read in the grass database was
> not correct : always the first one. Because in fact, it didn't manage to
> delete the temporary file before writing the new one..
>
> The best I can do, is give you part of the code as as you can see the
> context.
>
> Thanks for the code improvements with file.remove and unlink !
>
> Cheers,
>
> Mathieu
>
> for (id in 1:length(grids)) {
>
> (...)
>
> ## SMOOTH
>   print(">> SMOOTH")
>   smoothed = smooth(sieved, 9)
>   values(smoothed)[values(smoothed) > 0] = 1
>   values(smoothed)[values(smoothed)==0] = NA
>
>   ## THIN
>   # WRITE TO GRASS
>   writeRAST(as(smoothed, "SpatialGridDataFrame"), "bdtopo", zcol="layer",
> flags=c("overwrite"))
>   execGRASS("g.region", rast="bdtopo")
>   execGRASS("r.thin", parameters=list(input="bdtopo", output="thin"),
> flags=c("overwrite"))
>   thin = raster(readRAST("thin"))
>
>   ## CLEAN
>   execGRASS("r.to.vect", parameters=list(input="thin", output="thin",
> feature="line"), flags=c("overwrite"))
>   execGRASS("v.clean", parameters=list(input="thin", tool="rmdangle",
> output="clean", thresh=20), flags=c("overwrite"))
>
>   # REMOVE FILES
>   G = gmeta()
>   tmpDir = file.path(G$GISDBASE, G$LOCATION_NAME, G$MAPSET, ".tmp")
>   file.remove( list.files(path = tmpDir, Pattern ="^clean.*$", full.names
> = TRUE))
>
>   # READ TO SPDF
>   execGRASS("v.out.ogr", parameters=list(input="clean", type="line",
> dsn=tmpDir, olayer="clean", format="ESRI_Shapefile"))
>   cleaned = readOGR(tmpDir, "clean")
>
>   ## EXPORT
>   writeOGR(cleaned, "OUT/CLEANED", id, "ESRI Shapefile")
> }
>
> 2015-10-27 14:11 GMT+01:00 Rainer M Krug :
>
>> Mathieu Rajerison  writes:
>>
>> > Dear Mr Krug,
>> >
>>
>> I'll Cc the R-sig-geo list in to keep the conversation there for future
>> reference and for other users who might have the same problem.
>>
>> > Sorry in advance for my english.
>>
>> Don't worry abut your English - it is perfectly fine.
>>
>> >
>> > My question was not about debugging, but about a potentiel generic
>> > behaviour of spgrass6.
>> >
>> > I should have reformulated

Re: [R-sig-Geo] spgrass6, readVECT and overwriting ?

[Mathieu Rajerison]

Mr Krug,


It is difficult to reproduce an example from the spearfish dataset as it
doesn't natively contain data in it close to the one I used.

To give you th context, my goal is to detect all blue features : rivers
from a topographic map, ie 3-band RGB raster.

With RStoolbox, I perform a Spectral Angle Classification of my data using
sample points located on (near)blue pixels

With GRASS, I then thin, vectorize and clean the result.

As the process must be performed on a very large scale, I split it using
grids where sample RGB parts of the topo map are extracted with
gdal_translate. Everything is executed within a loop. At the end of the
loop, I export the resulted cleant shapefile in a folder. At the end, I
have as many shapefiles as processing grids.

That's why I use temporary vector datasets inside GRASS that are
overwritten at each loop instead of storing each individual result. I
noticed that in my process, the vector data read in the grass database was
not correct : always the first one. Because in fact, it didn't manage to
delete the temporary file before writing the new one..

The best I can do, is give you part of the code as as you can see the
context.

Thanks for the code improvements with file.remove and unlink !

Cheers,

Mathieu

for (id in 1:length(grids)) {

(...)

## SMOOTH
  print(">> SMOOTH")
  smoothed = smooth(sieved, 9)
  values(smoothed)[values(smoothed) > 0] = 1
  values(smoothed)[values(smoothed)==0] = NA

  ## THIN
  # WRITE TO GRASS
  writeRAST(as(smoothed, "SpatialGridDataFrame"), "bdtopo", zcol="layer",
flags=c("overwrite"))
  execGRASS("g.region", rast="bdtopo")
  execGRASS("r.thin", parameters=list(input="bdtopo", output="thin"),
flags=c("overwrite"))
  thin = raster(readRAST("thin"))

  ## CLEAN
  execGRASS("r.to.vect", parameters=list(input="thin", output="thin",
feature="line"), flags=c("overwrite"))
  execGRASS("v.clean", parameters=list(input="thin", tool="rmdangle",
output="clean", thresh=20), flags=c("overwrite"))

  # REMOVE FILES
  G = gmeta()
  tmpDir = file.path(G$GISDBASE, G$LOCATION_NAME, G$MAPSET, ".tmp")
  file.remove( list.files(path = tmpDir, Pattern ="^clean.*$", full.names =
TRUE))

  # READ TO SPDF
  execGRASS("v.out.ogr", parameters=list(input="clean", type="line",
dsn=tmpDir, olayer="clean", format="ESRI_Shapefile"))
  cleaned = readOGR(tmpDir, "clean")

  ## EXPORT
  writeOGR(cleaned, "OUT/CLEANED", id, "ESRI Shapefile")
}

2015-10-27 14:11 GMT+01:00 Rainer M Krug :

> Mathieu Rajerison  writes:
>
> > Dear Mr Krug,
> >
>
> I'll Cc the R-sig-geo list in to keep the conversation there for future
> reference and for other users who might have the same problem.
>
> > Sorry in advance for my english.
>
> Don't worry abut your English - it is perfectly fine.
>
> >
> > My question was not about debugging, but about a potentiel generic
> > behaviour of spgrass6.
> >
> > I should have reformulated my question in a simpler manner : does
> spgrass6
> > allow overwriting an existing dataset when using readVECT function ?
>
> The functions readVect / RAST are using temporary files - you are right,
> but they should delete these upon exit. So if they do nbot do this, it
> is a bug.
>
> >
> > 1) no version information about R and GRASS (there are several GRASS 6
> > versions)
> > I've just put it at the end of this post.
> > 2) no session information i.e. package versions used
> > the latest spgrass6
> > 3) no information about the OS (I guess it is windows based on the
> > Source path above)
> > you're right
>
> OK - thanks.
>
> > 4) no reproducible example which could be easily done by using the GRASS
>
> As I said - if you manually have to delete the temporary files, then
> there is a bug in the package. Could you please provide a reproducible
> example, using one of the example data sets from GRASS, so that I can
> look at it?
>
> > 5) as far as I can see, in the newest version of spgrass6 there is no
> > function readVECT()
> > Yes, there is, as you will see on page 14 of the documentation.
>
> Sorry - overlooked it.
>
> >
> > Finally, I solved the problem by using the following hint, "clean" being
> > the name of my temporary shapefile dataset :
> >   l = list.files("D:/GRASSDB/paca/mapset/.tmp", "^clean.*$"); l =
> > file.path("D:/GRASSDB/paca/mapset/.tmp", l)
> >   lapply(l,  file.remove)
> >   execGRASS

[R-sig-Geo] spgrass6, readVECT and overwriting ?

[Mathieu Rajerison]

Hi,

I'm using spgrass6 and I use readVECT function in a loop.

Using it causes an error because it doesn't overwrite the the temporary
shapefile..

> thin = readVECT(vname="zoneBDTOPOthin")
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERROR :Layer  already exists in OGR data source
'D:/GRASSDB/paca/mapset/.tmp'
OGR data source with driver: ESRI Shapefile
Source: "D:/GRASSDB/paca/mapset/.tmp", layer: "zoneBDTO"
with 770 features
It has 3 fields
Warning message:
running command 'v.out.ogr.exe -e -c input=zoneBDTOPOthin2 type=line
layer=1 dsn=D:/GRASSDB/paca/mapset/.tmp olayer=zoneBDTO
format=ESRI_Shapefile' had status 1


I wonder how to specify to spgrass that I want to overwrite the file in
readVECT, although I know using execGRASS and v.out.ogr with the
appropriate overwrite flags would do the job.

Thanks in advance for the answer..

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Complete traversals: Postman problem with shape files

[Mathieu Rajerison]

Hi,

Maybe you could have a look at v.net.salesman from GRASS.

Please note that GRASS can be reached while using R with spgrass06 package.

Best,

mathieu

2015-03-28 19:01 GMT+01:00 Mike :

> Hi all,
>
> I'm relatively new to the listserv, but learning a lot.  Please excuse me
> if this is a semi-related question.
>
> As part of a public health project to survey an area for tobacco
> retailers, we're attempting to investigate a systematic, complete, and
> reasonably efficient (not necessarily perfect) traversal of all the roads
> in a certain catchment area.  This is an application of the classic
> "postman" problem - traversing all the edges of a graph - related to the
> canonical Traveling Postman Problems / TSPs.   Specifically, we'd be taking
> regions of VA and breaking that up into volunteer catchment areas - a
> county or two - to drive the area.
> http://en.wikipedia.org/wiki/Route_inspection_problem
>
> I've a proof of concept plot and code below that takes a shape file of
> regions (in this case, of VA counties) and traverses them in a reasonably
> efficient way... but this uses centroids of shapes and creates a distance
> matrix off of that, instead of working with the road shape file itself,
> which is currently stumping me.
>
> (Below is me muddling through a test case using the counties as a whole,
> code and picture.)
>
> I know this is not at all a trivial problem.  Any guidance or expertise
> here?  I could make accessible more useful test cases/datasets (like a
> single county's roads) if that were helpful.  Or perhaps there are packages
> that do this out of the box besides what I'm using?
>
> Best wishes and much respect,
> Mike Fliss
> UNC Epidemiology PhD student.
>
>
> [image: Inline image 1]
>
> library("maptools")
> library("spdep")
> library("sp")
> library("TSP")
> #data("USCA312")
>
> VAcounties = readShapeSpatial("E:/Mike/GIS/Base shape files/VA Counties/VA
> Counties.shp")
> plot(VAcounties, border="grey")
> neighbors <-poly2nb(VAcounties)
> plot(neighbors, coordinates(VAcounties), line="grey", add=TRUE)
>
> VApoints = read.csv("E:/Dropbox/Community/Driving VA/vapoints.csv")
> VA.SPDF = SpatialPointsDataFrame(VApoints[,2:3], VApoints[1], proj4string
> = CRS("+proj=longlat"))
>
> vamatrix <- read.csv("E:/Dropbox/Community/Driving VA/vacountymatrix.csv",
> stringsAsFactors=FALSE)
> va2 = as.matrix(vamatrix)
> dimnames(va2) = list(names(vamatrix), names(vamatrix))
>
> vatsp = TSP(va2)
> vatsp=insert_dummy(vatsp, label="cut")
>
> #tour = solve_TSP(vatsp, method="repetitive_nn", control=list(start=1))
> tour = solve_TSP(vatsp, method="2-opt", control=list(rep=100))
> path = cut_tour(tour, "cut")
>
> head(labels(path))
>
> plot(VAcounties, border="grey")
> plot(VA.SPDF, pch=16, cex=.4, col="red", add=T)
> path_line = SpatialLines(list(Lines(list(Line(VA.SPDF[path,])), ID="1")))
> plot(path_line, add=T, col="black")
> points(VA.SPDF[c(head(path,1), tail(path,1)),], pch = 19, col = "black")
>
> --
> ---
> Mike Dolan Fliss, MSW
> mike.dolan.fl...@gmail.com
> UNC-CH Epidemiology PhD student
> NC public health advocate!
> ---
> "We work on ourselves in order to help others,
> but also we help others in order to work on ourselves."
>   - Pema Chodron
>
> “Upon this gifted age, in its dark hour
> Rains from the sky a meteoric shower
> Of facts…they lie, unquestioned, uncombined.
> Wisdom enough to leach us of our ill
> Is daily spun; but there exists no loom
> To weave it into a fabric.”
>   - Edna St. Vincent Millay, 1939
>
> There are two kinds of people in the world:
> Those who can extrapolate from incomplete data.
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] points sampled along a Line don't seem to intersect it

[Mathieu Rajerison]

Ok thanks for the answer.

Your explanations are very technical as they go deep into the package's
code, and GEOS specifications.

Finally, I used spatstat::project2Segment to get my points on the lines.

Mathieu

2015-03-23 20:30 GMT+01:00 Roger Bivand :

> On Mon, 23 Mar 2015, Barry Rowlingson wrote:
>
>  On Mon, Mar 23, 2015 at 12:07 PM, Mathieu Rajerison
>>  wrote:
>>
>>> pts <- spsample(rds.un, type="regular", n=100)
>>>
>>
>> Your points are at least 3.9x10^-17 units from the lines:
>>
>> gDistance(pts, rds.un)
>> [1] 3.908727e-17
>>
>> I can only think this is due to a difference in how spsample and how
>> rgeos interpolate lines from points. FAQ 7.31 in disguise?
>>
>> I suspect you could use gBuffer on your points with a buffer width of
>> 1e-10 to get a reliable intersection, but that would be a line segment
>> (or segments if the point is close to a junction) but you could
>> possibly just take a single point from that line geometry and in the
>> worst case scenario you're only 1e-10 out...
>>
>>
> Had the fuzzyMM package been written in a modularised way (it places
> vehicle GPS positions on OSM road lines, and needs velocity etc.), it might
> have given a starting point. As Barry says, finding which road is closest
> to each point is trivial (I projected too, to get more digits), but finding
> the nearest GEOS point on the road isn't, because GEOS discretises using
> getScale() to go to an integer grid. spsample() is also very pushed to get
> n here, so the # pts is way under.
>
> If anyone would like a brain teaser, this should be attractive!
>
> Roger
>
>
>>
>> Barry
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
> e-mail: roger.biv...@nhh.no
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] points sampled along a Line don't seem to intersect it

[Mathieu Rajerison]

Hi,


I have angle values affected to each segment of a road network and I'd like
to get the angles for sampling points.

The problem is that my points don't seem to intersect with my roads network.

Here is the code and the data is attached to this mail.

> library(rgeos)

> pts <- spsample(rds.un, type="regular", n=100)

> gIntersects(pts, rds.un)
[1] FALSE


Any idea ?

Best,

Mathieu


rds.un.rda
Description: Binary data
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] customized focal mean

[Mathieu Rajerison]

Hi,

Maybe you could look at raster::focal function and w (matrix of weights
arguments)
http://www.inside-r.org/packages/cran/raster/docs/focal

"w
matrix of weights (the moving window), e.g. a 3 by 3 matrix with values 1;
see Details. The matrix does not need to be square, but the sides must be
odd numbers. If you need even sides, you can add a column or row with
weights of zero"

2015-03-12 6:37 GMT+01:00 hariom :

> Hello list,
>
> I want to customized focal mean in which my concern is to replace mean
> value
> with the matrix such as:
>
> *2 2 2* 2
> *2 1 2* 1
> *1 2 4* 2
> 2 1 2 2
>
> first generate the mean of 3*3 matrix and then replace with entire matrix
> and it will work throughout the raster data.
> first matrix mean is 2 and  it will replace and continue.
>
>  2 *2 2 2*
>  2 *2 2 1*
>  2 *2 2 2*
>
>
> Please help me i will appreciate for your time and concern
>
> Hariom singh
> Research scholar IIT Roorkee
>
>
>
>
> --
> View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/customized-focal-mean-tp7587890.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Polygons VS MultiLineString

[Mathieu Rajerison]

Hello,

My advice is : why don't you reclass your raster in categories according to
raster values intervals, then you polygonize your reclassified raster ?

It would be simpler and you would have a topologically correct layer.

That's the way I do it with GRASS but you all the functions also are
included in the raster package.

Mat



2015-03-09 12:44 GMT+01:00 Antonio Rodriges :

> The solution I've devised so far
>
> contours <- rasterToContour(to_c, levels = c(18, 33))
>
> res <- lapply(slot(contours, "lines"), function(x) lapply(slot(x,
> "Lines"), function(y) slot(y, "coords")))
> polies <- lapply(res, function (x) lapply(x, function (y) Polygon(y)))
> pp_18 <- Polygons(polies[[1]], ID = "18")
> pp_33 <- Polygons(polies[[2]], ID = "33")
> all_pp<- SpatialPolygons(c(pp_18, pp_33))
> df <- data.frame(value = c(18, 33), row.names = c("18", "33"))
> SPDF <- SpatialPolygonsDataFrame(all_pp, df)
>
> writeOGR(SPDF, jsfile, layer="", driver="GeoJSON")
>
> This is in case you have only two values 18 and 33. If you have more,
> create list and create Polygons also using lapply or for loop
>
> 2015-03-05 11:49 GMT+03:00 Antonio Rodriges :
> > Mike,
> >
> > Thank you, that was a great example; I reproduced it.
> >
> > However, if I am satisfied with the rasterToContour output, I think
> > there should be a way to directly convert it polygons. I am seeking
> > this method.
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Generate N points along a SpatialLines object

[Mathieu Rajerison]

Ok nice !

I did'nt even notice spsample could work on SpatialLines object !..

2015-02-26 14:40 GMT+01:00 Edzer Pebesma :

> Mathieu, sp::spsample samples SpatialLines objects, and has an argument
> type = "regular".
>
> On 02/26/2015 02:24 PM, Mathieu Rajerison wrote:
> > Hi List,
> >
> >
> > I try to generate N points along a a continental coastline with
> > spatstat::pointsOnLines
> >
> > If I give a small number of points (like N = 10), then 931 points are
> > returned. That's because 931 segments are generated from my initial
> > coastline data.
> >
> > I'd like to know how to generate exactly N points, equally spaced, from a
> > SpatialLines object, without GRASS, nor QGIS
> >
> >
> > Thanks in advance for any answer,
> >
> > Mathieu
> >
> >   [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> 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/
> Spatial Statistics Society http://www.spatialstatistics.info
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Generate N points along a SpatialLines object

[Mathieu Rajerison]

Hi List,


I try to generate N points along a a continental coastline with
spatstat::pointsOnLines

If I give a small number of points (like N = 10), then 931 points are
returned. That's because 931 segments are generated from my initial
coastline data.

I'd like to know how to generate exactly N points, equally spaced, from a
SpatialLines object, without GRASS, nor QGIS


Thanks in advance for any answer,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to segment or split a spatial line in R

[Mathieu Rajerison]

Hi,

With spatstat, you can find the mid point of line segments
http://lojze.lugos.si/~darja/software/r/library/spatstat/html/midpoints.psp.html

But I don't know any package that would split your line into equal segments.

One method would be to arrange equally spaced points along your line with
spatstat::pointsOnLinesthen reconstruct lines from these points.

There may be a better method than mine...

2014-10-23 16:10 GMT+02:00 Manuel Spínola :

> Dear list members,
>
> How to segment or split a spatial line in shorter equal segments, and also,
> how to get the mid point of ecah segment.
>
> Best,
>
> Manuel
>
> --
> *Manuel Spínola, Ph.D.*
> Instituto Internacional en Conservación y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspin...@una.ac.cr
> mspinol...@gmail.com
> Teléfono: (506) 2277-3598
> Fax: (506) 2237-7036
> Personal website: Lobito de río <
> https://sites.google.com/site/lobitoderio/>
> Institutional website: ICOMVIS 
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Kernel Density Estimation with Border Bias Correction

[Mathieu Rajerison]

Hello,

Some people here might be interested by this :
https://github.com/ripleyCorr/Kernel_density_ripley

"This page proposes some R codes to compute the kernel density estimates of
two-dimensional data points, using an extension of Ripley's circumference
method to correct for border bias"

Best regards,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Anselin Local Moran's I with R

[Mathieu Rajerison]

Hello,

See a loot at spdep, especially localmoran function
http://www.inside-r.org/packages/cran/spdep/docs/localmoran

Best,

Mathieu

2014-09-09 4:32 GMT+02:00 David Romero :

> Hello,
>
> Could somebody help me with the method to compute the Cluster and Outlier
> Analysis using spdep and obtain Local I Index, Z-scores and cluster type
> like with the Arcgis tool.
> My data are points with temperature value.
>
> Thank you,
>
> David Romero
> phD student
> Instituto de Geografía
> UNAM
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to do an anamorphosis on R?

[Mathieu Rajerison]

Hi,

There are many solutions, the most popular of which is using the ScapeToad
software outside of R

See http://stackoverflow.com/questions/9319597/cartogram-choropleth-map-in-r
and
http://spatial.ly/2013/06/r_activity/

Mat


2014-02-24 16:24 GMT+01:00 Gilles Benjamin Leduc :

> Hi all,
>
> I would like to make an anamorphosis map, id est, a map where some other
> distance measurement are used instead geographical distance. In my case I
> would like to plot an anamorphosis of Iceland using genetical distence
> insteand of geographical distence.
> I have got .shp maps and equivalents.
> How to make such a distortion?
> Best regards
> Benjamin
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Shortest-path (railway) distance calculation using R

[Mathieu Rajerison]

Hi,

You could use igraph package to do this
See Barry Rowlingson's method  : http://rpubs.com/geospacedman/routing

Otherwrise, with OSM data, you can use osmar package.
But I don't think you can coerce a SpatialLines object to an osmar one.
Maybe if you convert your data into OSM XML data using GDAL first ?
http://osmar.r-forge.r-project.org/

Hope this helps,

Mathieu


2013/10/24 Daisuke Murakami 

> Thank you very much!
> Daisuke
>
> > Have you looked at http://cran.r-project.org/web/views/Spatial.html,
> > specifically the `rgeos` package?
> >
> > Cheers,
> > Roman
> >
> >
> >
> >
> > On Wed, Oct 23, 2013 at 11:03 AM, Daisuke Murakami <
> > murak...@sk.tsukuba.ac.jp> wrote:
> >
> >> Hi everyone,
> >>
> >> I try to calculate railway distances from a railway station to the other
> >> stations on R.
> >> I have shp files of the stations and railway.
> >>
> >> Does anyone know useful R package/function for the distance calculation?
> >>
> >> Sincerely,
> >> Daisuke
> >>
> >> ___
> >> R-sig-Geo mailing list
> >> R-sig-Geo@r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> >
> >
> > --
> > In God we trust, all others bring data.
> >
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] voronoi (thiessen) polygons inside irregular area

[Mathieu Rajerison]

Hi,

spatstat::dirichlet can do that
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/spatstat/html/dirichlet.html

Note that to get them as SpatialPolygons thereafter, you've got to do this :
spatstat.options("gpclib")
= TRUE

best

mathieu


2013/10/2 Milan Cisty 

> I would like to draw voronoi or thiessen polygons inside my watershed
> polygon which I have in shape file. I found only functions for drawing
> thiessen polygons in rectangular area - please could you let me know, if
> function for finding thiessen polygons exists also for irregular
> boundaries?
> I will need also what will be the area of each polygon.
>
> Thanks
>
> Milan
>
>
>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Scaling SpatialPolygons

[Mathieu Rajerison]

Hi,

Also, see maptools::elide which gives some functions for that.

You could take some material from this git if you want:
https://gist.github.com/datagistips/5262410#file-cartograms-r

with 2 functions: one using buffer, another one using elide

best,

mathieu


2013/7/17 Rizzo, Michael 

> Robert,
> Thank you!  That did exactly what I needed.
>
> Mike Rizzo
>
> -Original Message-
> From: Robert J. Hijmans [mailto:r.hijm...@gmail.com]
> Sent: Wednesday, July 17, 2013 12:01 PM
> To: Rizzo, Michael
> Cc: r-sig-geo@r-project.org
> Subject: Re: [R-sig-Geo] Scaling SpatialPolygons
>
> Mike, You can create a small inside (negative) buffer using rgeos and use
> that. Robert
>
> On Wed, Jul 17, 2013 at 9:17 AM, Rizzo, Michael 
> wrote:
> > I have a group of counties in a SpatialPolygonsDataFrame and would like
> to know if a function exists to scale them.  Many of the counties share a
> common border, and I want to color each county's border differently to
> distinguish them by individual categories.  The problem is the border
> colors overlap for adjacent counties.   I thought if I could scale the
> individual polygons down to 95% of their original size, the boundaries of
> adjacent counties would appear side by side without the overlap.
> >
> > I created my own function to do this but concave polygons don't scale
> 100% correctly and still overlap each other a little so I am hoping there
> may be a pre-defined function that would scale the county polygons
> correctly.
> >
> > Thank you,
> > Mike Rizzo
> >
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] what's new in asdar latest edition?

[Mathieu Rajerison]

ok, edzer, thanks for the reply!

congratulations for this book that let me discover the wonderful world of
spatial stats


2013/7/2 Roger Bivand 

> On Mon, 1 Jul 2013, Edzer Pebesma wrote:
>
>  Mathieu,
>>
>> The second edition of "Applied Spatial Data Analysis with R" (asdar2)
>> seems to be shipping now, and contains the following changes to chapters:
>> - ch 6, "Classes for spatio-temporal Data", replaces the old ch 6,
>> "Customising spatial data classes and methods". The old ch 6 is now a
>> vignette in sp:
>> http://cran.r-project.org/web/**packages/sp/vignettes/csdacm.**pdf<http://cran.r-project.org/web/packages/sp/vignettes/csdacm.pdf>
>> - ch 7, "Spatial Point Pattern Analysis", pays, compared to the first
>> edition, more attention to package spatstat
>> - ch 9, Modelling Areal Data is a merge of the old ch 9, "Areal data
>> and spatial correlation" and ch 10, "Modelling areal data"
>> - Chapter 10 (former 11) has been expanded to include recent
>> developments in Bayesian modelling, including BayesX and R-INLA
>> packages, and new examples. Spatio-temporal models for disease mapping
>> are described now as well.
>>
>
> The book website, http://www.asdar-book.org, now has the second edition
> code, code+data bundles, and errata. We'd be grateful for reports on
> mistakes and things that don't work to add to the errata, so that
> users/readers can keep running smoothly. We don't yet have the code for
> INLA, WinBUGS or GRASS up yet - to follow.
>
> Roger
>
>
>
>> The second edition contains full colour figures, many new examples, and
>> addresses several new packages, in particular rgeos and spacetime.
>> References were updated or added. Code examples have been updated, e.g.
>> the use of overlay() was removed altogether in favour of over(),
>> aggregate() has been explained, and spatial selections based on
>> intersection are now done by x[y,] with y a Spatial* object.
>>
>> After the spatial workshop at UseR!, next week (Jul 9, 17:45 CEST), a
>> little celebration will take place to launch the second edition --
>> everyone is invited! Roger and Virgilio will be there in person.
>> https://twitter.com/**SpringerStats/statuses/**350673280354615298<https://twitter.com/SpringerStats/statuses/350673280354615298>
>>
>> Book site:
>> http://www.springer.com/**statistics/life+sciences%2C+**
>> medicine+%26+health/book/978-**1-4614-7617-7<http://www.springer.com/statistics/life+sciences%2C+medicine+%26+health/book/978-1-4614-7617-7>
>>
>> The first edition of the book has now appeared in a chinese translation:
>> http://www.openbookdata.com.**cn/Book/MzA2NDAyOF9PU0dB.html<http://www.openbookdata.com.cn/Book/MzA2NDAyOF9PU0dB.html>
>>
>> Best regards,
>>
>> On 05/02/2013 12:17 PM, Mathieu Rajerison wrote:
>>
>>> Hello,
>>>
>>>
>>> I've already bought asdar's previous edition.
>>>
>>> It's an excellent book that I recommend. As a beginner in spatial stats,
>>> I
>>> learnt a lot and the content seemed accessible for me.
>>>
>>> Apparently, there's a new edition of asdar book, which I don't know a lot
>>> about.
>>>
>>> I'd like to know some reasons I'd have to buy the newest edition. New
>>> chapters, content?
>>>
>>>
>>> Mathieu
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> __**_
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>
>>>
>>
>>
> --
> Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: roger.biv...@nhh.no
>
>
> __**_
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spatial networks

[Mathieu Rajerison]

Hi,

Good news.

For me, what was missing for the moment with spatial tasks was something
that could coerce linear spatial objects to (i)graph objects so as to pull
out all the potentiel of igraph package.

I can remember some attempts by Barry Rownlingson with a buildTopo function
http://rpubs.com/geospacedman/routing

What I find constraining in the SpatialLinesNetwork function is that it
only accepts singleLine SpatialLines objects, althoug the fact is that the
spatial data we'll use will mainly be composed of roads, rivers, etc...that
aren't, natively, single lines

Maybe we could suggest a function that "explodes" multi-line SpatialLines
objects while preserving the attributes (since the attributes can drive
some graph analysises)

Also, one more thing I think of is snapping: let's suppose I have two
points representing buildings and I'd like to compute the shortest path
between them. With spatstat, we can get the edge connecting the node (the
building) to the network, so we're not far from executing the graph stats.

To sum up, I'd sugget:
- accepting multi-line SpatialLines objects
- building a new network given extra/outer nodes

For the moment, to do the basic graph stats, I'd still use GRASS . I would
be happy to do all this stuff with R. And I think we're really not that far
from that..



2013/6/25 Forrest Stevens 

> Hi Edzer, this is definitely very cool.  This may make some recent
> work I've been doing on spatial and social network analyses of
> academic co-authorship easier (early code is here if you're
> interested:  http://refnet.r-forge.r-project.org/ ).  Thanks for
> providing this to the community!
>
> Sincerely,
> Forrest
>
> On Sun, Jun 23, 2013 at 10:20 PM, Thomas Adams  wrote:
> > Edzer,
> >
> > That's really pretty cool stuff; I think I can find some applications for
> > it!
> >
> > Cheers!
> > Tom
> >
> >
> > On Sun, Jun 23, 2013 at 5:57 PM, Edzer Pebesma <
> > edzer.pebe...@uni-muenster.de> wrote:
> >
> >> I did some simple things trying to combine sp and igraph objects to
> >> create spatial networks, and published this at:
> >>
> >> http://rpubs.com/edzer/spatialnetworks
> >>
> >> I'd be happy to share the markdown document with you; it seems to me
> >> that the rpubs don't allow you to simply download & run all the commands
> >> in the html.
> >>
> >> Anyway, comments & discussion welcome, as usual.
> >> --
> >> Edzer Pebesma
> >> Institute for Geoinformatics (ifgi), University of Münster
> >> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> >> 83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
> >>
> >> ___
> >> R-sig-Geo mailing list
> >> R-sig-Geo@r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> > [[alternative HTML version deleted]]
> >
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
>
> --
> Forrest R. Stevens
> Ph.D. Candidate, QSE3 IGERT Fellow
> Department of Geography
> Land Use and Environmental Change Institute
> University of Florida
> www.clas.ufl.edu/users/forrest
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] SpatialLines to Graph object - computing the longest path - misconstruction

[Mathieu Rajerison]

Hello,


I have a lines shapefile that I converted to a graph object, following the
method in http://www.mail-archive.com/r-sig-geo@r-project.org/msg05836.html

I have determined the farthest nodes of the graph.

Then, I tried to determine the diameter's path  but the result is odd as
the attached file shows.

The problem may come from the SpatialLines to graph conversion..But I can't
see where the problem is...

Any help would be greatly appreciated...

I've put the example data I used if you want to give it a try

Here is the code:

# READ
l <- readOGR("IN/testline.shp", "testline")
l.break <- gIntersection(l, l)
l.split <- SpatialLines(lapply(1:length(l.break@lines[[1]]@Lines),
function(x) Lines(list(l.break@lines[[1]]@Lines[[x]]), ID=as.character(x

#GRAPH
edges =
do.call(rbind,lapply(1:length(l.split),function(x){as.vector(t(l.split@lines
[[x]]@Lines[[1]]@coords))}))
lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)

froms = paste(edges[,1],edges[,2])
tos = paste(edges[,3],edges[,4])
g = graph.edgelist(cbind(froms,tos),directed=FALSE)
E(g)$weight=lengths

# DIAMETER
sel <- get.diameter(g)
diam <- l.split[sel-1]
plot(diam)


Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Skeletonization Process - Eliminating Adjacent Segments

[Mathieu Rajerison]

Hi,

Finally, I inspired myself from a code by Barry Rownlingson intially
dedicated to calculate the shortest path.
But in my cas, it consisted in calculating the longest one.

Here is the final code with an image as attached file. In the image, you
see the initial shape, the voronoi lines in rainbow colors and the
simplified median Line in bold black.

library(rgdal)
library(spatstat)
library(rgeos)
library(maptools)
library(igraph)

setwd("F:/METHODES/1305_SKELETON/")

pol <- readOGR("IN/polygone.shp", "polygone")

# FUNCTION

findMedianLine <- function(pol, eps, tol) {

  l <- as(as(pol, "SpatialLines"), "psp")

  pts <- pointsOnLines(l, eps=eps)

  vor <- dirichlet(pts)


  # LINES EXTRACTION

  vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")

  vor.l <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines

  vor.split <- SpatialLines(lapply(1:length(vor.l), function(x)
Lines(list(vor.l[[x]]), ID=as.character(x

  vor.l.in <- vor.l.split[gContains(pol, vor.split, byid=T), ]


  # GRAPH CONSTRUCTION

  edges = do.call(rbind,lapply(vor.l.in@lines
,function(ls){as.vector(t(ls@Lines[[1]]@coords))}))
  lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)

  froms = paste(edges[,1],edges[,2])
  tos = paste(edges[,3],edges[,4])

  g = graph.edgelist(cbind(froms,tos),directed=FALSE)
  E(g)$weight=lengths


  # LONGEST PATH

  sel <- get.diameter(g)

  vor.sel <- vor.l.in@lines[sel-1]

  spine <- SpatialLines(vor.sel)

  spine.spl <- gSimplify(gLineMerge(spine), tol=tol)

  return(spine.spl)
}


myMedianLine <- findMedianLine(pol, eps=1000, tol=1000)


2013/5/17 Mathieu Rajerison 

> Hi,
>
>
> I've written a code to skeletonize an area shape.
>
> The example could be a watershed from which we'd like to compute the
> "guessed" stream, or find the median axis of a shape.
>
> I can find my median line but there are some adjacent segments connected
> to it. The length criteria is not appropriate to deal with these adjacent
> segments as some tiny ones could be part of to the median line.
>
> Has anyone idea to deal with that?
>
> I apoplogize in advance if it doesn't work in all cases (holes, rings,
> etc...).
>
> library(rgdal)
> library(spatstat)
> library(rgeos)
> library(maptools)
>
> setwd("F:/METHODES/1305_SKELETON/")
>
> pol <- readOGR("IN/polygone.shp", "polygone")
>
> # FUNCTION
>
> findMedianLine <- function(pol, tol) {
>
>   l <- as(as(pol, "SpatialLines"), "psp")
>
>   pts <- pointsOnLines(l, eps=1)
>
>   vor <- dirichlet(pts)
>
>   vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")
>
>   vor.lines <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines
>
>   vor.l.split <- SpatialLines(lapply(1:length(vor.lines), function(x)
> Lines(list(vor.lines[[x]]), ID=as.character(x
>
>   medianLine <- gLineMerge(vor.l.split[gContains(pol, vor.l.split,
> byid=T), ])
>
>   medianLine.spl <- gSimplify(medianLine, tol=tol); plot(medianLine.spl)
>
>   return(medianLine.spl)
> }
>
>
> myMedianLine <- findMedianLine(pol, tol=1000)
>
> plot(pol); plot(myMedianLine, add=T)
>
>
>
>
>
<>___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Skeletonization Process - Eliminating Adjacent Segments

[Mathieu Rajerison]

Hi,


I've written a code to skeletonize an area shape.

The example could be a watershed from which we'd like to compute the
"guessed" stream, or find the median axis of a shape.

I can find my median line but there are some adjacent segments connected to
it. The length criteria is not appropriate to deal with these adjacent
segments as some tiny ones could be part of to the median line.

Has anyone idea to deal with that?

I apoplogize in advance if it doesn't work in all cases (holes, rings,
etc...).

library(rgdal)
library(spatstat)
library(rgeos)
library(maptools)

setwd("F:/METHODES/1305_SKELETON/")

pol <- readOGR("IN/polygone.shp", "polygone")

# FUNCTION

findMedianLine <- function(pol, tol) {

  l <- as(as(pol, "SpatialLines"), "psp")

  pts <- pointsOnLines(l, eps=1)

  vor <- dirichlet(pts)

  vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")

  vor.lines <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines

  vor.l.split <- SpatialLines(lapply(1:length(vor.lines), function(x)
Lines(list(vor.lines[[x]]), ID=as.character(x

  medianLine <- gLineMerge(vor.l.split[gContains(pol, vor.l.split, byid=T),
])

  medianLine.spl <- gSimplify(medianLine, tol=tol); plot(medianLine.spl)

  return(medianLine.spl)
}


myMedianLine <- findMedianLine(pol, tol=1000)

plot(pol); plot(myMedianLine, add=T)

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Finding the Circumscribing Circle of a polygon

[Mathieu Rajerison]

Thanks Michael for the answer.

How come that I didn't find this solution on google?...When you use rseek
and type circumscribing circle, you don't find any mention to the tripack
package..Maybe it's because the fnction has been developed lately? the
latest version of tripack being from february 2013

Anyway, it's nice to have a solution, now.

What's the best source when you look for a function to accomplish something?


2013/5/16 Michael Sumner 

> Just FYI, there is  tripack::circumcircle, here's an example with data
> in the maptools package.
>
> library(maptools)
>
> data(wrld_simpl)
>
> poly <- wrld_simpl[wrld_simpl$NAME == "Australia", ]
>
> ## use coercion to get raw vertices (long winded, but avoids @ usage)
> coords <- coordinates(as(as(poly, "SpatialLines"), "SpatialPoints"))
>
>
> library(tripack)
> cc<- circumcircle(coords[,1], coords[,2], plot = TRUE)
> plot(poly, add = TRUE, col = "grey")
>
>
> On Fri, May 17, 2013 at 7:39 AM, JLong  wrote:
> > Hey Mathieu and Greg,
> >
> > Great ideas, and thanks for the code Mathieu.
> > I can't imagine a polygon shape where the 2nd step would ever by
> necessary?
> > Can you elaborate on why that might be.
> >
> > FYI, Mathieu I'm attempting to compute the same shape/linearity metric as
> > you.
> >
> > Cheers,
> > Jed
> >
> >
> >
> >
> > -
> > Jed Long
> > PhD Candidate
> > SPAR Lab, Dept. of Geography
> > University of Victoria, Canada
> > --
> > View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/Finding-the-Circumscribing-Circle-of-a-polygon-tp7582372p7583598.html
> > Sent from the R-sig-geo mailing list archive at Nabble.com.
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> Michael Sumner
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Modis tiles with different resolutions

[Mathieu Rajerison]

Hi,

Try raster::resample then raster::mosaic


2013/5/9 lancelot 

> Dear all,
>
> I am processing large sets of MODIS data that come to me in sinusoidal
> projection and .lan (ERDAS 7.5) format. I have to mosaic different tiles
> and here comes the problem I meet: resolutions are not the same.
>
> Basically I'm looping through these steps:
>
> library(raster)
> library(sp)
> wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
> sinu <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181
> +units=m +no_defs"
> ## MA (h22v10) tile extent in longlat
> e <- extent(c(40.6187, 53.2072, -20, -10))
> ## transform to sinusoidal
> e <- as.data.frame(t(bbox(e)))
> coordinates(e) <- ~ s1 + s2
> projection(e) <- wgs84
> e <- spTransform(e, CRS(sinu))
> e <- coordinates(e)
> e <- extent(c(e[,1], e[,2]))
> ## First raster
> r <- raster("c:/Temp/RtmpyQ3t4R/**MA050107.lan")
> projection(r) <- sinu
> extent(r) <- e
>
> at the end of the process:
>
> > ## tile MA (h22v10)
> > MAdlst[[1]]
> class   : RasterLayer
> dimensions  : 1200, 1200, 144  (nrow, ncol, ncell)
> resolution  : 1318.567, 926.6254  (x, y)
> extent  : 4244214, 5826494, -2223901, -951  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
> +b=6371007.181 +units=m +no_defs
> data source : c:\Temp\RtmpyQ3t4R\MA050107.**lan
> names   : MA050107
> values  : -32768, 32767  (min, max)
>
> > ## tile MB (h22v11)
> > MBdlst[[1]]
> class   : RasterLayer
> dimensions  : 1200, 1200, 144  (nrow, ncol, ncell)
> resolution  : 1611.032, 926.6254  (x, y)
> extent  : 4099193, 6032431, -3335852, -2223901  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
> +b=6371007.181 +units=m +no_defs
> data source : c:\Temp\RtmpyQ3t4R\MB050107.**lan
> names   : MB050107
> values  : -32768, 32767  (min, max)
>
>
> Am I doing something wrong, and how can I mosaic these rasters with
> different resolutions?
>
> All the best,
>
> Renaud
>
>
> > sessionInfo()
> R version 3.0.0 (2013-04-03)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
> LC_MONETARY=French_France.1252 LC_NUMERIC=C
> [5] LC_TIME=French_France.1252
>
> attached base packages:
> [1] grid  stats graphics  grDevices utils datasets  methods
> base
>
> other attached packages:
> [1] rgdal_0.8-9 maptools_0.8-23 foreign_0.8-53  lattice_0.20-15
> mapdata_2.2-2   maps_2.3-2  raster_2.1-25
> [8] sp_1.0-9
>
> loaded via a namespace (and not attached):
> [1] compiler_3.0.0 fortunes_1.5-1 tools_3.0.0
>
>
>
>
>
>
>
>
> Le 08/05/2013 21:45, Francisco Zambrano a écrit :
>
>> Robert,
>>
>> Sorry, now I realize that calc works in chunk if is needed. And I imagine
>> that calc is more efficient than do it the chunks by myself.
>>
>> Thanks,
>>
>> Francisco
>>
>>
>> 2013/5/8 Robert J. Hijmans
>>
>>  Francisco,
>>>
>>> I think what you want to do can be accomplished with 'calc' (as was
>>> pointed out in another thread).
>>>
>>> library(raster)
>>> in<- stack(list.files(pattern=*..**tif$))
>>> out<- calc(in, fun=function(x) lowess(x,f=7/n,iter=10)$y)
>>>
>>>
>>> But to answer your question, I think you could do:
>>>
>>> in<- stack(list.files(pattern=*..**tif$))
>>> # and, here is the trick
>>> out<- brick(in, values=FALSE)
>>>
>>> Also, in your loop you could do:
>>>
>>> getValues(in, ...
>>> writeValues(out, 
>>>
>>>
>>> Best, Robert
>>>
>>> On Tue, May 7, 2013 at 7:35 AM, Francisco Zambrano
>>> wrote:
>>>
 Hi everybody,

 I'm working with MODIS data MOD13Q1 and the NDVI index between 2000-

>>> today.
>>>
 I have been trying to do a function, but I have a problem when stack the
 310 files and then I tried to use the writeStart function. I had an
 error
 because the startFuction work with rasterLayer or RasterBrick. But, when

>>> I
>>>
 tried to make the brick with the rasterStack (>brick(rasterStack)) the
 operation took too much time (11 hours)

 My question, is what is the most efficiently way to make a rasterBrick

>>> from
>>>
 multiple raster (geoTIFF) files?

 save the rasterStack and then open it with brick?

 My code is something like:

  out<-stack(list.files(pattern=***..tif$))
> out<-brick(mod) #too much time (11 hours appox)
>

>out<-writeStart(out,filename,**overwrite=TRUE)
>bs<-blockSize(out)
>pb<-pbCreate(bs$n)
>for (i in 1:bs$n){
  v<-getValues(out,row=bs$row[i]**,nrows=bs$nrows[i])
  s<-t(sapply(apply(v,1,lowess,**f=7/n,iter=10),"[[","y"))
  out<-writeValues(out,s,bs$row[**i])
 pbStep(pb,i)
}

 Kind Regards

 Francisco Zambrano Bigiarini

 from Chile, latin america

  [[alternative HTML version deleted]]

 __**_
 R-sig-Geo mailing 

Re: [R-sig-Geo] Determining the shared boundaries between polygons/polyshapes

[Mathieu Rajerison]

Hi,

Maybe, using rgeos and intersecting two copies of your vegetation layer?
Kind of gIntersection(vegLayer, vergLayer, byid=TRUE) ?


2013/5/5 Alastair Potts 

> Good day,
>
> I would like to determine the percentage/length of shared boundaries
> between a set of polygons (which are, in this case, mapped vegetation
> units).
>
> This has been dealt with in the r-sig-geo list before (
> https://stat.ethz.ch/pipermail/r-sig-geo/2010-August/008917.html);
> however,
> the solution involved porting across to GRASS (which is something I would
> like to avoid given my beginner status of GIS in R).
>
> I was wondering if this problem is now solvable in R, and if so, where
> should I look?
>
> Thanks in advance,
>
> Regards,
> Alastair
>
> --
> Dr. Alastair J. Potts
> Claude Leon Postdoctoral Fellow
> Botany Department
> Nelson Mandela Metropolitan University
>
> Cell #: 082 491-7275
> Office #: 041 504-4375 (Mon, Thurs, Fri)
> Fax #: 086 273-2675
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] what's new in asdar latest edition?

[Mathieu Rajerison]

Hello,


I've already bought asdar's previous edition.

It's an excellent book that I recommend. As a beginner in spatial stats, I
learnt a lot and the content seemed accessible for me.

Apparently, there's a new edition of asdar book, which I don't know a lot
about.

I'd like to know some reasons I'd have to buy the newest edition. New
chapters, content?


Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Distance converted in koordinates

[Mathieu Rajerison]

Hi,

Maybe you could write this to a file, then read this file using space
delimiter or TAB?


2013/4/17 Endri Raco 

> Hi group,
> I am stucked in a problem like this.
> I have distances from a bus station in km generated in R language like
>
>
> 0.  0.02725136  1.07634471  1.15963225  1.71421571  2.54945626
> 4.29135102  4.53532958  4.58512710  4.86466833  5.24266630  5.63505465
>
> I need to convert this distances in coordinates in order to reflect
> them as points in map
>
> Any ideas about this?
>
> Please help
>
> Regards
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] drawing ellipses to display uncertainty around estimates of animal locations

[Mathieu Rajerison]

Hi,

sp:bubble function?


2013/4/17 Mónica Cordeiro de Almeida e Silva 

> Dear All,
>
> I'm using switching state-space models (SSSM) to analyze Argos-derived
> locations from whales tracked with satellite tags. I'd like to draw
> ellipses to display the uncertainty at each location estimated by the SSSM,
> using only 3 lat/long pairs of values provided by the model: the mean
> location and the locations corresponding to the 2.5 and 97.5 percentiles
> (ie. the 0.95 probability interval). Each ellipse should be centered at the
> mean value and its limits should be defined by the coordinates of the 2.5
> and 97.5 percentiles.
>
> Any hints on how to do this in R? Any suggestions would be welcome!
>
> Cheers,
>
> Mónica
>
>
> [[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster::distance, "raster has a longitude/latitude CRS, but coordinates do not match that"

[Mathieu Rajerison]

Hi,

I've got a vector layer of the world that I took from Natural Earth.

I tried to get a distance raster from the coastline. With a raster
resolution of 10, no problem, I get the distance raster. But if I decrease
the resolution to 5, I get the error:

"1: In .couldBeLonLat(x) :
  raster has a longitude/latitude CRS, but coordinates do not match that
2: In raster:::.couldBeLonLat(x) :
  raster has a longitude/latitude CRS, but coordinates do not match that"


Here is the code:

library(rgdal)
library(raster)

world <- readOGR("F:/DATAS/NATURAL
EARTH/ne_110m_admin_0_countries_lakes/ne_110m_admin_0_countries_lakes.shp",
"ne_110m_admin_0_countries_lakes")

r <- raster(extent(world)); res(r) <- 5; proj4string(r) <-
CRS("+init=EPSG:4326")
world.r <- rasterize(world, r)
dist <- distance(world.r)


Here is my session info:

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C   LC_TIME=French_France.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] raster_2.0-12 rgdal_0.7-18  sp_0.9-99

loaded via a namespace (and not attached):
[1] grid_2.15.1lattice_0.20-6 tools_2.15.1


Thanks for any help,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] randomly sample points

[Mathieu Rajerison]

Sorry, I meant polygonize via raster::rasterToPolygons

2013/3/8 Mathieu Rajerison 

> Hi,
>
> You could rasterize you layer then apply spsample function on each cell.
>
>
> 2013/3/8 Ross Ahmed 
>
>> The following code randomly samples cells within a raster:
>>
>> myRaster <- raster(ncol=10, nrow=10)
>> values(myRaster) <- rnorm(ncell(myRaster))
>> plot(myRaster)
>> library(dismo)
>> xy <- randomPoints(myRaster, n=100)
>> points(xy)
>>
>> Is it possible to randomly sample points within each cell?
>>
>> Thanks
>> Ross
>>
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] randomly sample points

[Mathieu Rajerison]

Hi,

You could rasterize you layer then apply spsample function on each cell.

2013/3/8 Ross Ahmed 

> The following code randomly samples cells within a raster:
>
> myRaster <- raster(ncol=10, nrow=10)
> values(myRaster) <- rnorm(ncell(myRaster))
> plot(myRaster)
> library(dismo)
> xy <- randomPoints(myRaster, n=100)
> points(xy)
>
> Is it possible to randomly sample points within each cell?
>
> Thanks
> Ross
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Spatial points to lines/polygons

[Mathieu Rajerison]

Hi,

You could create the spatialpolygons from scratch using this example:
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/sp/html/SpatialPolygons-class.html

I noticed you tried to create an alphahull object. Maybe you should have a
look at this code published by Barry Rowlingson:
http://rpubs.com/geospacedman/alphasimple

Also, you've got the rgeos readWKT function to create geometries from Well
Known Text Geometry coordinates.

Best,

Mathieu


2013/2/11 Gustaf Granath 

> Hi,
> I have mapped structures with a GPS (walked along their edges) and Im now
> trying to make polygons from my coordinates. The GPS data are a bit fuzzy
> and I want to establish the edge of the structure with a smoother or
> similar. However, I have not been able to find an easy way to do this in R.
> This seems odd to me and I might have missed something. I have so far tried
> doing lines with spline (xspline) or polygons with the alphahull package,
> but they are tricky to work with (usually overfit).   Any pointers are
> welcome.
>
> Gustaf
>
>
> Example code, testing the alphahull package:
> require(alphahull)
> set.seed(1)
> xx <-c(1:10,8,5,2)
> xx <- xx + rnorm(39,0,0.4)
> yy <-c(1:5,5:1,1,1,1)
> yy <-yy + rnorm(39,0,0.4)
> plot( xx ~ yy)
> xy <-data.frame("x" = xx, "y" = yy)
> coords <- coordinates(xy)
>
> # alpha-shape
> x.as <- ashape(coords, alpha=5)
>
> # alpha-hull
> x.ah <- ahull(coords, alpha=5)
>
> #plot
> plot(x.as)
> plot(x.ah)
>
> --
> Gustaf Granath (PhD)
> Post doc
> McMaster University
> Email: gustaf.gran...@gmail.com
>
> __**_
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Finding the Circumscribing Circle of a polygon

[Mathieu Rajerison]

Hello,


Here is a code based on Greg's second idea (1st step: find the furthest
points of the polygon, 2nd step: find the circumscribed circle of a
triangle)

I wrote a function that computes the smallest circumscribing circle of a
triangle, and one that searches for the furthest vertices of a polygon.

The second function could help you find the median axis of a shape.

For the smallest circumscribing, I just took the formula from wikipedia.

You'll find a dummy example illustrating the search for the circle of a
triangle (2nd step).

You'll also find a spatial example including the search of the furthest two
points (1st step).

I've tested this code on a few countries from natural earth. It worked
quite well. Except sometimes, the polygon overlaps the circle.

Also, it performs the operation on the first entity. So, if the shape is of
type multipolygon, the algorithm will only work on the first polygon it
finds.

I think I now have a sufficient material to calculate the linearity index
but don't hesitate if you see something wrong in my code or have advice for
some algorithm improvements.


#
## DUMMY EXAMPLE ##
#

# http://en.wikipedia.org/wiki/Circumscribed_circle

findCircumTriangle <- function(tricoords) {
vert.x <- tricoords[,1]
vert.y <- tricoords[,2]

D <- 2*(vert.x[1]*(vert.y[2]-vert.y[3])
  + vert.x[2]*(vert.y[3]-vert.y[1])
  + vert.x[3]*(vert.y[1]-vert.y[2]))

circum.x <- ((vert.x[1]^2 + vert.y[1]^2)*(vert.y[2]-vert.y[3])
 +(vert.x[2]^2 + vert.y[2]^2)*(vert.y[3]-vert.y[1])
 +(vert.x[3]^2 + vert.y[3]^2)*(vert.y[1]-vert.y[2])) / D

circum.y <- ((vert.x[1]^2 + vert.y[1]^2)*(vert.x[3]-vert.x[2])
 +(vert.x[2]^2 + vert.y[2]^2)*(vert.x[1]-vert.x[3])
 +(vert.x[3]^2 + vert.y[3]^2)*(vert.x[2]-vert.x[1])) / D

dx <- (vert.x[1] - circum.x)
dy <- (vert.y[1] - circum.y)
r <- sqrt(dx^2 + dy^2)

return(list(x=circum.x, y=circum.y, radius=r))
}


coords <- cbind(c(1,0,3), c(.5,2,3))
circum <- findCircumTriangle(coords); circum

plot( 0:5, 0:5, type='n', xlab='', ylab='', asp=1)
polygon(coords)
theta <- seq(0,2*pi, length.out=250)
polygon( circum$x + circum$r*cos(theta),
 circum$y + circum$r*sin(theta), border='blue')


##
## SPATIAL EXAMPLE ##
##

library(rgdal)

countries <- readOGR("D:/DATAS/NATURAL
EARTH/ne_110m_admin_0_countries_lakes/ne_110m_admin_0_countries_lakes.shp",
"ne_110m_admin_0_countries_lakes")

country <- countries[24, ]
coords <- country@polygons[[1]]@Polygons[[1]]@coords

findMidBetweenFurthest <-  function(polycoords) {
  d <-  as.matrix(dist(polycoords))
  dmax <- max(d)

  first <- which(d==dmax, arr.ind=T)[1,1]
  second <- which(d==dmax, arr.ind=T)[1,2]
  mid <- apply(cbind(polycoords[first, ], polycoords[second, ]), 1,
FUN=mean)

  return(list(first=first , second=second, coords=mid, radius=abs(dmax/2)))

}

# 1st Step
mid <- findMidBetweenFurthest(coords)

# PLOTTING
theta <- seq(0,2*pi, length.out=250)
polygon( mid$coords[1] + mid$radius*cos(theta),
 mid$coords[2] + mid$radius*sin(theta), border='blue' )

d <- tail(as.matrix(dist(rbind(coords, mid$coords))), 1)
furthest <- which(d == max(d), arr.ind=T)[1,1]

# 2nd STEP
if (max(d) > mid$radius) {
  coordsTriangle <- rbind(coords[furthest, ], coords[mid$first, ],
coords[mid$second, ])
  circum <- findCircumTriangle(coordsTriangle)
} else {
  circum <- list(x=mid$coords[1], y=mid$coords[2], radius=mid$radius)
}

# PLOTTING
plot(country)
points(rbind(coords[mid$first, ], mid$coords, coords[mid$second, ]),
col="red")
theta <- seq(0,2*pi, length.out=250)
polygon( circum$x + circum$r*cos(theta),
 circum$y + circum$r*sin(theta), border='blue')


2013/1/26 Greg Snow <538...@gmail.com>

> The "dists" variable represents the distance between the points and the
> circle measured along the radius of the circle.  Positive distances mean
> that the point/vertex is inside the circle and negative distances mean that
> the point is outside of the circle.  Without some type of penalty the
> optimal circle would generally have some of the points inside and some
> outside (this would be a regression like model fitting the best circle
> instead of the best line).  Since we want all the points/vertices to be in
> the circle we add the penalty for any points outside the circle.  This is
> the penalty that worked for my examples, some other penalty may work
> better.
>
> For a more analytic solution see my second post.  That idea could possibly
> use some improvement as well, it is just what I thought of between waking
> up and actually getting out of bed (then refined somewhat while waiting for
> train into work to come).  It makes sense to me that it should work for
> most cases, maybe even all cases, but I have not worked out a formal proof
> of that.
>
>
> On Sat, Jan 26, 2013 at 3:57 AM, Berry Boessenkool <
> berryboessenk...@hotmail.com> wrote:
>
> >
> > Hi Greg,
> >
> > why did you add in partic

Re: [R-sig-Geo] Finding the Circumscribing Circle of a polygon

[Mathieu Rajerison]

Thanks a lot for your answer.

I wasn't familiar with optimization methods.

One question I asked myself was how to determine the optimal center
coordinates and radius to start the optimization with.

So, based on your code:
- I used the centroid as the center
- half the size of the biggest side of the bounding box as radius

I got one problem. After optimizing, the circle still touched the shape.
So, I made a loop that increased the radius until the shape didn't touch
the circle.

I simply translated your method to make it work on SpatialPolygons object.

I used a sample from natural earth countries data as a sample to test your
method and it worked quite well.

Thanks again for your answer.

The code, for the moment based on Greg Snow's solution (any advice
appreciated)

library(rgeos)

myfun <- function(theta) {
  x <- theta[1]
  y <- theta[2]
  r <- theta[3]

  dists <- r - sqrt( (vert.x-x)^2 + (vert.y-y)^2 )
  sum( dists^2 ) + 10*sum(dists<0)
}

spP = countries[13, ] # sampling a country from Natural Earth Data
spP = SpatialPolygons(list(Srs1))

# POLYGON COORDINATES
vert.x = spP@polygons[[1]]@Polygons[[1]]@coords[, 1]
vert.y = spP@polygons[[1]]@Polygons[[1]]@coords[, 2]

# OPTIMIZATION
coords <- coordinates(spP)
ini.r <- min(diff(t(bbox(spP
out <- optim(c(coords[1], coords[2], ini.r), myfun)

# CIRCLE
theta <- seq(0,2*pi, length.out=250)
coordPol = cbind(out$par[1] + out$par[3]*cos(theta), out$par[2] +
out$par[3]*sin(theta))
coordPol = rbind(coordPol, coordPol[1, ])
circle = SpatialPolygons(list(Polygons(list(Polygon(coordPol)), "pol")))

# BUFFER CIRCLE UNTIL NO INTERSECTION
i = 0
inc = 0.01
d = out$par[3] * inc
while(!gContains(circle, spP)) {
  d = d * (1 + inc)
  circle <- gBuffer(circle, width = radius)
  i = i+1
  print(i)
}

# PLOTTING
plot(spP)
plot(circle, add=T)

# LINEARITY INDEX
lin = 1 - (gArea(spP)/gArea(circle))


2013/1/25 Greg Snow <538...@gmail.com>

> Here is one possibility to get you close:
>
> vert.x <- c(0,0,1,1)
> vert.y <- c(0,1,1,0)
>
> myfun <- function(theta) {
> x <- theta[1]
>  y <- theta[2]
> r <- theta[3]
>
> dists <- r - sqrt( (vert.x-x)^2 + (vert.y-y)^2 )
>  sum( dists^2 ) + 10*sum(dists<0)
> }
>
> out <- optim( c(0,0,5), myfun )
> out
>
>
> plot( -1:2, -1:2, type='n', xlab='', ylab='', asp=1 )
> polygon(vert.x, vert.y)
>
> theta <- seq(0,2*pi, length.out=250)
> polygon( out$par[1] + out$par[3]*cos(theta),
> out$par[2] + out$par[3]*sin(theta), border='blue' )
>
> tmp.x <- rnorm(25)
> tmp.y <- rnorm(25)
> tmp <- chull(tmp.x,tmp.y)
>
> vert.x <- tmp.x[tmp]
> vert.y <- tmp.y[tmp]
>
> plot(vert.x, vert.y, xlim=c(-3,3), ylim=c(-3,3))
> polygon(vert.x,vert.y)
>
> out <- optim( c(0,0,3), myfun )
> out
>
> theta <- seq(0,2*pi, length.out=250)
> polygon( out$par[1] + out$par[3]*cos(theta),
> out$par[2] + out$par[3]*sin(theta), border='blue' )
>
>
>
> On Thu, Jan 24, 2013 at 7:19 AM, Mathieu Rajerison <
> mathieu.rajeri...@gmail.com> wrote:
>
>> I made a little mistake. The formula is 1 - area(shape) / area(circum),
>> not 1 - area(circum) / area(shape)
>>
>> 2013/1/24 Mathieu Rajerison 
>>
>> > Hi,
>> >
>> >
>> > I know spatstat::incircle to calculate the inscribed circle of a shape.
>> >
>> > I'd like to know if there was a function to calculate the circumscribing
>> > circle of a polygon.
>> >
>> > The area of this circumscribing circle could then be used to calculate a
>> > linearity index: 1 - area(circum) / area(shape)
>> >
>> > For a triangle, the circumcenter is at the intersection of the different
>> > perpendicular bisectors of its sides and each corner would be at the
>> same
>> > distance to this center, distance being known as circumradius.
>> >
>> > For a polygon, it seems to me it's a lot more complicated.
>> >
>> > It's more a geometrical problem than a spatial one but if anyone has an
>> > advice, I'd be glad of it..
>> >
>> >
>> > Mathieu
>> >
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
> --
> Gregory (Greg) L. Snow Ph.D.
> 538...@gmail.com
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Finding the Circumscribing Circle of a polygon

[Mathieu Rajerison]

I made a little mistake. The formula is 1 - area(shape) / area(circum),
not 1 - area(circum) / area(shape)

2013/1/24 Mathieu Rajerison 

> Hi,
>
>
> I know spatstat::incircle to calculate the inscribed circle of a shape.
>
> I'd like to know if there was a function to calculate the circumscribing
> circle of a polygon.
>
> The area of this circumscribing circle could then be used to calculate a
> linearity index: 1 - area(circum) / area(shape)
>
> For a triangle, the circumcenter is at the intersection of the different
> perpendicular bisectors of its sides and each corner would be at the same
> distance to this center, distance being known as circumradius.
>
> For a polygon, it seems to me it's a lot more complicated.
>
> It's more a geometrical problem than a spatial one but if anyone has an
> advice, I'd be glad of it..
>
>
> Mathieu
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Finding the Circumscribing Circle of a polygon

[Mathieu Rajerison]

Hi,


I know spatstat::incircle to calculate the inscribed circle of a shape.

I'd like to know if there was a function to calculate the circumscribing
circle of a polygon.

The area of this circumscribing circle could then be used to calculate a
linearity index: 1 - area(circum) / area(shape)

For a triangle, the circumcenter is at the intersection of the different
perpendicular bisectors of its sides and each corner would be at the same
distance to this center, distance being known as circumradius.

For a polygon, it seems to me it's a lot more complicated.

It's more a geometrical problem than a spatial one but if anyone has an
advice, I'd be glad of it..


Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] breaking lines at intersection but without GRASS

[Mathieu Rajerison]

Lu,

Thanks for the answer. I had already looked at shp2graph which looked
really fine but that I had trouble installing it properly.

It's nice if we can install it easily now. I'll have to give it a try

To note, also, that the osmar package which provides access to osm data and
coercion to Spatial* Objects can also convert data to igraph objects.

Best,

Mathieu

2012/10/12 Barry Rowlingson 

> I did all that from scratch:
>
> http://www.rpubs.com/geospacedman/routing
>
> not sure how it compares with shp2graph... I'll have to try sometime!
>
> Barry
>
>
> On Fri, Oct 12, 2012 at 11:47 AM, binbin lu  wrote:
> > Hi,
> >
> > The package shp2graph has been repacked and the bug has been fixed. It
> will
> > be fine to install for the latest version of R (2.15). Hope it's still
> > useful for you.
> >
> > Cheers.
> >
> > Binbin
> >
> > 2012/6/27 Mathieu Rajerison 
> >
> >> Hi,
> >>
> >>
> >> My final goal is to generate a graph from a road network. I didn't
> manage
> >> to install shp2graph correctly, so I want to do it from scratch.
> >>
> >> Is there a way to break lines at intersection but without rgrass?
> >>
> >> Also, apart from shp2graph which seems a very good package, does one
> exist
> >> that would turn my road network into a graph?
> >>
> >>
> >> Thanks!
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> ___
> >> R-sig-Geo mailing list
> >> R-sig-Geo@r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Finding the highest order of neighbors

[Mathieu Rajerison]

Hi rgeo list!


Considering a neighbor object, how can I find the highest order of
neighbors and which node corresponds to it?

I think of using spdep::nblag iteratively until I find no change between
the lower order and the actual one, but I find the method a bit "dirty".

Would you have any idea?


Thanks in advance,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] simplify a polygon window

[Mathieu Rajerison]

Hi Sara,

You can't avoid this when simplifying.

But you could find the window that results from the intersection from your
simplified one and the original with intersect.owin.

This window will not overlap the original one.


2012/10/4 sara martino 

> Hi,
>
> I am working with a point pattern whose window is a spatial polygon which
> describe the cost of Norway.
>
> I need to get a smoother version of such polygon which wraps the original
> one.
> I have tried to use the function
>
> > simplify.owin(my.owin,8)
> but this does not work since some parts of the original polygon will lay
> outside the simplified one.
>
> Does anyone know how to do that?
> thanks a lot
> sara
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Problem aggregating municipalities

[Mathieu Rajerison]

Hi,

Maybe you could ckeck, first, which geometries are not valid, using
gIsValid and using par. byid=TRUE
v <- gIsValid(vM.spdf, byid=TRUE)

Then, you could fix the topology problems.


2012/9/30 Brigitte Hogan 

> I am having a problem with unionSpatialPolygons() in
> library(maptools). I am working with a shapefile of Mexico
> municipalities found here:
> http://mapserver.inegi.org.mx/data/mgm/
> under Marco Geoestadístico Municipal 2009 Versión 4.1,  Áreas
> Geoestadísticas Municipales (38.8 Mb).
>
> I need to combine the "Federal District" municipalities into one
> polygon. I am a new user of the sp and rgeos packages, running R
> version 2.13 under Windows 7. When I try to combine the districts, I
> get an error similar to one in the archive:
> https://stat.ethz.ch/pipermail/r-sig-geo/2011-August/012445.html
>
> I tried to run checkPolygonsHoles() as suggested, but I cannot figure
> out how to convert the SpatialPolygonsDataFrame or SpatialPolygons
> object into a Polygons object.  The solution to the error in the
> archive was to reset the significant digits. However, when I tried, I
> still ended up with the same error message. An excerpt from my code is
> below.
>
>
> ## Original Code
> library(maptools)
> library(rgdal)
> library(rgeos)
> file3 <- "C:/Desktop/MUNICIPIOS"
> vM.spdf <- readShapePoly(file3, IDvar=NULL)
> class(vM.spdf)
># SpatialPolygonsDataFrame
> vM.spdf@data[which(vM.spdf@data$CVE_ENT=="09"),]
> v.id <- vM.spdf@data$CVE_ENT=="09"   # IDs
> which rows to merge
> newM <- unionSpatialPolygons(vM.sp, IDs=v.id)  # merge
>
> # Returns the error:
> # Error in TopologyFunc(groupID(spgeom[ids[[i]]], id[ids[[i]]]),
> names(ids)[i],  :
> # TopologyException: found non-noded intersection between LINESTRING
> (3.47895e+006 551354,
> # 3.479e+006 551323) and LINESTRING (3.479e+006 551323, 3.479e+006
> 551364) at 3.479e+006 551323
>
> ## Trying checkPolygonsHoles()
> checkPolygonsHoles(vM.spdf)#
> Error ... not a Polygons object
> vM.sp <- SpatialPolygons(vM.spdf@polygons, vM.spdf@plotOrder,
> proj4string=vM.spdf@proj4string)
> vM.p   <- Polygons(vM.sp)
> # Error in as.list.default(X) : no method for coercing this S4 class
> to a vector
> vM.p  <- Polygons(list(vM.sp))
> vM.p <- Polygons(vM.sp@polygons) # Error
> in Polygons(list(vM.sp)) : srl not a list of Polygon objects
>
> ## Trying to reset sig difs
> getScale()
>   # returns '1e+08'
> setScale(1e+09)
> newM <- unionSpatialPolygons(vM.sp, IDs=v.id)  # same
> error
>
>
>
> Thanks for any help you can give me.
>
> Brigitte
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] NA values in rasterStack

[Mathieu Rajerison]

Hi Julio,

To put NA values on your raster for values equal to 3:
r[r==3] <- NA

2012/9/20 Julio Oliveira 

> Dear,
>
> I have two rasterStack (NDVI and B) with 232 layers each.
>
> NDVI
> class   : RasterBrick
> dimensions  : 10, 11, 110, 232  (nrow, ncol, ncell, nlayers)
> resolution  : 0.00225, 0.00225  (x, y)
> extent  : -48.345, -48.32025, -20.3625, -20.34  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
> values  : in memory
> min values  : 0.46 0.64 0.58 0.46 0.30 0.20 0.19 0.24 0.20 0.22 ...
> max values  : 0.91 0.95 0.94 0.94 0.88 0.76 0.70 0.82 0.77 0.78 ...
> layer names : MOD13Q1.A2000353.h13v11.005.2006346191943.250m_16_days_NDVI,
> MOD13Q1.A2001001.h13v11.005.2008270060743.250m_16_days_NDVI,
> MOD13Q1.A2001017.h13v11.005.2006358162247.250m_16_days_NDVI,
> MOD13Q1.A2001033.h13v11.005.2008271231754.250m_16_days_NDVI,
> MOD13Q1.A2001049.h13v11.005.2007002215512.250m_16_days_NDVI,
> MOD13Q1.A2001065.h13v11.005.2007007103033.250m_16_days_NDVI,
> MOD13Q1.A2001081.h13v11.005.2008278010133.250m_16_days_NDVI,
> MOD13Q1.A2001097.h13v11.005.2007018110448.250m_16_days_NDVI,
> MOD13Q1.A2001113.h13v11.005.2008284095248.250m_16_days_NDVI,
> MOD13Q1.A2001129.h13v11.005.2007027112124.250m_16_days_NDVI, .
>
>
> IMAGE B
> The image B has values 0,1 or  3 for each layer.
>
> class   : RasterBrick
> dimensions  : 10, 11, 110, 232  (nrow, ncol, ncell, nlayers)
> resolution  : 0.00225, 0.00225  (x, y)
> extent  : -48.345, -48.32025, -20.3625, -20.34  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
> values  : in memory
> min values  : 1 0 0 0 0 0 0 0 0 0 ...
> max values  : 3 0 0 0 0 3 0 0 0 1 ...
> layer names :
> MOD13Q1.A2000353.h13v11.005.2006346191943.250m_16_days_pixel_reliability,
> MOD13Q1.A2001001.h13v11.005.2008270060743.250m_16_days_pixel_reliability,
> MOD13Q1.A2001017.h13v11.005.2006358162247.250m_16_days_pixel_reliability,
> MOD13Q1.A2001033.h13v11.005.2008271231754.250m_16_days_pixel_reliability,
> MOD13Q1.A2001049.h13v11.005.2007002215512.250m_16_days_pixel_reliability,
> MOD13Q1.A2001065.h13v11.005.2007007103033.250m_16_days_pixel_reliability,
> MOD13Q1.A2001081.h13v11.005.2008278010133.250m_16_days_pixel_reliability,
> MOD13Q1.A2001097.h13v11.005.2007018110448.250m_16_days_pixel_reliability,
> MOD13Q1.A2001113.h13v11.005.2008284095248.250m_16_days_pixel_reliability,
> MOD13Q1.A2001129.h13v11.005.2007027112124.250m_16_days_pixel_reliability
>
>
> First Question:
>
> I'll like select cells (ex.: values == 3) from difference layer in B  and
> put NA for each cells / layer in NDVI image.
>
> I tried the syntax:
>
> NDVI <- stackSelect(NDVI, !(B==3), recycle=TRUE, type = 'index')
> Error: nlx > nly is not TRUE
>
> Second Question:
> How put NA values for all the pixels with value == 3 on image B ?
>
> Thanks
>
> Julio
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] project2segment -> connecting SpatialLines: some connect, some not

[Mathieu Rajerison]

Hi,


This post could have followed the one about shortest path but I preferred
to make a separate post to make things clearer.

Given a set of points and a set of lines, I tried to make a graph of it.
The operation consisted, firstly, in getting the lines that connect each
point to the nearest point on the segment pattern

For this, I used spatstat::project2segment. spatstat is really a wonderful
package. I think it can be useful for snapping line operations.

I'm stuck because I manage to create some lines between the points and the
segments but some lines connect to the network while others don't.

When they don't, the distance between the end connecting segment point and
the network is extremely small. I put an image of it as attached file.

Here is the code, for those who'd like to give it a try and, eventually,
give an issue. Maybe I did it in a wrong way. It's based on Barry's initial
code.

library(spatstat)
library(rgeos)
library(maptools)

makeNest <- function(n) {
### generate a bunch of SpatialLines on (0,1) square
xy = cbind(runif(n), runif(n), runif(n), runif(n))
lines = Lines(apply(xy, 1, function(z) {
Line(rbind(z[1:2], z[3:4]))
}), ID = 1)
SpatialLines(list(lines))
}

connection2Lines <- function(lines, pts){

  sourcePts <- as(pts,"ppp")
  projectedPts <- project2segment(sourcePts, as(lines, "psp"))$Xproj
  bb <- union.owin(bounding.box(sourcePts), bounding.box(projectedPts))
  cnnLines <- as(psp(sourcePts$x, sourcePts$y, projectedPts$x,
projectedPts$y, window=bb), "SpatialLines")
  cnnLines <- spChFIDs(cnnLines, sapply(cnnLines@lines, function(x)
paste("connect", x@ID, sep="_")))
  proj4string(cnnLines) <- proj4string(lines)

  return(cnnLines)
}

n <- 5
pts <- SpatialPoints(cbind(runif(n), runif(n)))

gg = makeNest(10)

cnnLines <- connection2Lines(gg, pts)

plot(gg)
plot(cnnLines, col="red", add=T)

any(gIntersects(cnnLines, gg, byid=TRUE)==FALSE)
<>___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

Edzer,

I'm sorry.

The only control I did was checking if the resulting network looked right.
At a global extent, it was, but zooming in revealed some errors.

The reason why I didn't compute the shortest path calculation on my script
is that I didn't have the minimum knowledge in igraph to do it. That's why
I said I didn't test the script completely. I should have mentionned I
expected some feedback on the shortest path calculation.

Your feedback and the code published by Barry helped me do some deeper
inspection.

The fact is that the lines that must go from one point to the nearest point
on the lines network are well created but some connect to the network, some
not. For those who don't, the distance between the end segment point and
the line to which it is supposed is extremely small.

My intent wasn't to pollute this thread, rather to contribute to it, not as
an "expert", which I am not, but as a beginner. I'll be more cautious in
the future, with all the respect I owe to the R-sig-geo contributors who
participate on the list.

Mathieu

2012/9/11 Edzer Pebesma 

>
>
> On 09/10/2012 12:47 PM, Mathieu Rajerison wrote:
> > Barry,
> >
> > I looked at your code and I wondered if it would not be more convenient
> to
> > find the closest point on a an edge rather than the closest vertex?
> >
> > In this case, what technique would be the best? I thought of
> > spatstat::project2Segment that I used in my script (but Edzer says it
> > doesn't work - I have to check that)
>
> I only made a comment on sending untested scripts, I made no comment on
> whether project2Segment did what you want it to do.
>
> >
> >
> > Mathieu
> >
> >
> > 2012/9/10 Mathieu Rajerison 
> >
> >> It looks really nice.
> >>
> >> Thanks Barry!
> >>
> >>
> >> 2012/9/9 Barry Rowlingson 
> >>
> >>> I've just pushed a cleanup and a bunch of examples to RPubs:
> >>>
> >>> http://rpubs.com/geospacedman/routing
> >>>
> >>>  It uses the latest igraph!
> >>>
> >>>  It's all there, building the network and routing and plotting.
> >>>
> >>> Barry
> >>>
> >>> ___
> >>> R-sig-Geo mailing list
> >>> R-sig-Geo@r-project.org
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>
> >>
> >>
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
> http://www.52north.org/geostatistics  e.pebe...@wwu.de
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

Barry,

I looked at your code and I wondered if it would not be more convenient to
find the closest point on a an edge rather than the closest vertex?

In this case, what technique would be the best? I thought of
spatstat::project2Segment that I used in my script (but Edzer says it
doesn't work - I have to check that)


Mathieu


2012/9/10 Mathieu Rajerison 

> It looks really nice.
>
> Thanks Barry!
>
>
> 2012/9/9 Barry Rowlingson 
>
>> I've just pushed a cleanup and a bunch of examples to RPubs:
>>
>> http://rpubs.com/geospacedman/routing
>>
>>  It uses the latest igraph!
>>
>>  It's all there, building the network and routing and plotting.
>>
>> Barry
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

It looks really nice.

Thanks Barry!

2012/9/9 Barry Rowlingson 

> I've just pushed a cleanup and a bunch of examples to RPubs:
>
> http://rpubs.com/geospacedman/routing
>
>  It uses the latest igraph!
>
>  It's all there, building the network and routing and plotting.
>
> Barry
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

Hi,

My code isn't right as the Lines object has to be built from the beginning
by including connecting lines. Then the graph can be constructed.

I updated my initial thread on "Lines network to graph" with a new function
that corrects that:
http://r-sig-geo.2731867.n2.nabble.com/Lines-Network-to-Graph-td7580741.html


I haven't tested my function and hope it works well.

Some improvements can still be done in order to get an enriched function
for the conversion of a Lines network to a graph, which should be very
useful.

Mathieu

2012/9/1 Mathieu Rajerison 

> Hi,
>
> Following Barry's code and using spatstat::project2segment, I implemented
> the integration of extra nodes to the graph
>
> I tried to give different attributes to vertexes depending on whether they
> were on the lines or whether they're extra nodes, in order to facilitate
> queries but I don't find my code very clean.
>
> It's :
>
>   E(graph)$weight=lengths
>   nVtx <- length(V(graph))
>
>   # TYPE ATTRIBUTE
>   V(graph)$type <- "Line"
>   V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"
>
>
> Would you see any improvement to my code?
>
> Here is the "new" function based on Barry's one
>
> buildTopo <- function(lines, pts){
>   require(rgeos)
>   require(igraph)
>   require(spatstat)
>
>   # (1) LINEAR NETWORK
>   g = gIntersection(f,f)
>   edges1 = do.call(rbind,lapply(g@lines
> [[1]]@Lines,function(ls){as.vector(t(ls@coords))}))
>
>   # (2) POINT LOCATIONS
>   froms2 <- as(pts,"ppp")
>   tos2 <- project2segment(froms2, as(g, "psp"))$Xproj
>   edges2 <- cbind(froms2$x, froms2$y, tos2$x, tos2$y)
>
>   # MERGING BOTH ((1)+(2))
>   edges <- rbind(edges1, edges2)
>
>   # LENGTH CALCULATION
>   lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
>
>   # GRAPH
>   froms = paste(edges[,1],edges[,2])
>   tos = paste(edges[,3],edges[,4])
>
>   graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
>   E(graph)$weight=lengths
>   nVtx <- length(V(graph))
>
>   # TYPE ATTRIBUTE
>   V(graph)$type <- "Line"
>   V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"
>
>   return(graph)
> }
>
> m = readOGR(".","Line")
> pts <- readOGR(".", "Point")
> gg = buildTopo(m, pts)
> E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]
>
> I'm glad to have, now, a totally R topology function!
>
> Mathieu
>
> 2012/8/31 Mathieu Rajerison 
>
>> Wow, I'm really impressed by your geohacking skills. Thanks a lot.
>>
>> Yes, there is a GRASS::v.clean function that cuts lines.
>> Also, GRASS already includes allocation algorithms.
>>
>> We could also imagine performing shortest path calculations for extra
>> vertexes that aren't initally on the graph: shop locations, particular
>> events, houses...
>>
>> For this, spatstat::project2segment might be the appropriate function.
>>
>> Based on your initial code, I think we could write a function that would
>> combine a linear network with extra vertexes into a single graph. The
>> vertex labels could discriminate the intersection nodes from the thematic
>> ones.
>>
>>
>>
>> 2012/8/31 Barry Rowlingson 
>>
>>> On Fri, Aug 31, 2012 at 11:58 AM, Rowlingson, Barry
>>>  wrote:
>>>
>>> >  If I didn't have to be on a train in three hours I'd code this up...
>>>
>>>  Oh who am I kidding:
>>>
>>> buildTopo <- function(lines){
>>>   require(rgeos)
>>>   require(igraph)
>>>
>>>   g = gIntersection(lines,lines)
>>>   edges = do.call(rbind,lapply(g@lines
>>> [[1]]@Lines,function(ls){as.vector(t(ls@coords))}))
>>>   lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
>>>
>>>   froms = paste(edges[,1],edges[,2])
>>>   tos = paste(edges[,3],edges[,4])
>>>
>>>   graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
>>>   E(graph)$weight=lengths
>>>   return(graph)
>>>
>>> }
>>>
>>>  > m = readOGR(".","Line")
>>>  > gg = buildTopo(m)
>>>
>>> compute shortest distance between two vertices (note the -1 is because
>>> I'm still on old igraph with 0-indexing):
>>>
>>>  > E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]
>>>
>>> Edge sequence:
>>>
>>> [0]  0.81711193 47.57711634 -- -0.28546154 48.4287593
>>> [1]  2.2151139 46.49728033  -- 0.81711193 47.57711634
>>> [9]  1.57160852 45.36348513 -- 2.2151139 46.49728033
>>> [10] 1.39973834 45.06066625 -- 1.57160852 45.36348513
>>> [11] 1.29755246 44.88062446 -- 1.39973834 45.06066625
>>> [12] 0.91544563 43.67971729 -- 1.29755246 44.88062446
>>> [13] 0.88815229 43.43407719 -- 0.91544563 43.67971729
>>>
>>> I think this could be improved by adding the coordinates to the nodes.
>>> But anyway, 90% done.
>>>
>>> Barry
>>>
>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Lines Network to Graph

[Mathieu Rajerison]

Hi,

Finally, Barry came up with a solution.
http://r-sig-geo.2731867.n2.nabble.com/shortest-path-td7580832.html

I used his solution to take into account additional points, extra nodes.

For this, I connected points to their nearest point on the segment with
project2segment function. Then, I rebuilt the graph.

What I'd like to do next is assigning a diffferent attribute to nodes that
were extra vertexes using a column from the points dataset. They can be
shops, particular events, etc...

buildTopo <- function(lines, pts){
  require(rgeos)
  require(igraph)
  require(spatstat)

  # PROJECTION TO SEGMENT
  froms2 <- as(pts,"ppp")
  tos2 <- project2segment(froms2, as(lines, "psp"))$Xproj
  bb <- union.owin(bounding.box(froms2), bounding.box(tos2))
  connectLines <- as(psp(froms2$x, froms2$y, tos2$x, tos2$y, window=bb),
"SpatialLines")
  IDs <- sapply(connectLines@lines, function(x) paste("connect", x@ID,
sep="_"))
  connectLines <- spChFIDs(connectLines,IDs)
  proj4string(connectLines) <- proj4string(lines)

  # MERGING LINES
  Lines <- spRbind(lines, connectLines)

  # LINEAR NETWORK
  g = gIntersection(Lines,Lines)
  edges = do.call(rbind,lapply(g@lines
[[1]]@Lines,function(ls){as.vector(t(ls@coords))}))

  # LENGTH CALCULATION
  lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)

  # GRAPH
  froms = paste(edges[,1],edges[,2])
  tos = paste(edges[,3],edges[,4])

  graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
  E(graph)$weight=lengths

  return(graph)
}

Best,

Mathieu


2012/8/13 Mathieu Rajerison 

> Hi,
>
>
> Considering a linear network, how can I generate a graph of it that I
> could use in igraph? So that I could perform shortest paths and other
> calculations...
>
> The vertexes would be the intersections of lines.
>
>
> Thanks
>
> Mathieu
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

Hi,

Following Barry's code and using spatstat::project2segment, I implemented
the integration of extra nodes to the graph

I tried to give different attributes to vertexes depending on whether they
were on the lines or whether they're extra nodes, in order to facilitate
queries but I don't find my code very clean.

It's :

  E(graph)$weight=lengths
  nVtx <- length(V(graph))

  # TYPE ATTRIBUTE
  V(graph)$type <- "Line"
  V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"


Would you see any improvement to my code?

Here is the "new" function based on Barry's one

buildTopo <- function(lines, pts){
  require(rgeos)
  require(igraph)
  require(spatstat)

  # (1) LINEAR NETWORK
  g = gIntersection(f,f)
  edges1 = do.call(rbind,lapply(g@lines
[[1]]@Lines,function(ls){as.vector(t(ls@coords))}))

  # (2) POINT LOCATIONS
  froms2 <- as(pts,"ppp")
  tos2 <- project2segment(froms2, as(g, "psp"))$Xproj
  edges2 <- cbind(froms2$x, froms2$y, tos2$x, tos2$y)

  # MERGING BOTH ((1)+(2))
  edges <- rbind(edges1, edges2)

  # LENGTH CALCULATION
  lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)

  # GRAPH
  froms = paste(edges[,1],edges[,2])
  tos = paste(edges[,3],edges[,4])

  graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
  E(graph)$weight=lengths
  nVtx <- length(V(graph))

  # TYPE ATTRIBUTE
  V(graph)$type <- "Line"
  V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"

  return(graph)
}

m = readOGR(".","Line")
pts <- readOGR(".", "Point")
gg = buildTopo(m, pts)
E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]

I'm glad to have, now, a totally R topology function!

Mathieu

2012/8/31 Mathieu Rajerison 

> Wow, I'm really impressed by your geohacking skills. Thanks a lot.
>
> Yes, there is a GRASS::v.clean function that cuts lines.
> Also, GRASS already includes allocation algorithms.
>
> We could also imagine performing shortest path calculations for extra
> vertexes that aren't initally on the graph: shop locations, particular
> events, houses...
>
> For this, spatstat::project2segment might be the appropriate function.
>
> Based on your initial code, I think we could write a function that would
> combine a linear network with extra vertexes into a single graph. The
> vertex labels could discriminate the intersection nodes from the thematic
> ones.
>
>
>
> 2012/8/31 Barry Rowlingson 
>
>> On Fri, Aug 31, 2012 at 11:58 AM, Rowlingson, Barry
>>  wrote:
>>
>> >  If I didn't have to be on a train in three hours I'd code this up...
>>
>>  Oh who am I kidding:
>>
>> buildTopo <- function(lines){
>>   require(rgeos)
>>   require(igraph)
>>
>>   g = gIntersection(lines,lines)
>>   edges = do.call(rbind,lapply(g@lines
>> [[1]]@Lines,function(ls){as.vector(t(ls@coords))}))
>>   lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
>>
>>   froms = paste(edges[,1],edges[,2])
>>   tos = paste(edges[,3],edges[,4])
>>
>>   graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
>>   E(graph)$weight=lengths
>>   return(graph)
>>
>> }
>>
>>  > m = readOGR(".","Line")
>>  > gg = buildTopo(m)
>>
>> compute shortest distance between two vertices (note the -1 is because
>> I'm still on old igraph with 0-indexing):
>>
>>  > E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]
>>
>> Edge sequence:
>>
>> [0]  0.81711193 47.57711634 -- -0.28546154 48.4287593
>> [1]  2.2151139 46.49728033  -- 0.81711193 47.57711634
>> [9]  1.57160852 45.36348513 -- 2.2151139 46.49728033
>> [10] 1.39973834 45.06066625 -- 1.57160852 45.36348513
>> [11] 1.29755246 44.88062446 -- 1.39973834 45.06066625
>> [12] 0.91544563 43.67971729 -- 1.29755246 44.88062446
>> [13] 0.88815229 43.43407719 -- 0.91544563 43.67971729
>>
>> I think this could be improved by adding the coordinates to the nodes.
>> But anyway, 90% done.
>>
>> Barry
>>
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] shortest path

[Mathieu Rajerison]

Wow, I'm really impressed by your geohacking skills. Thanks a lot.

Yes, there is a GRASS::v.clean function that cuts lines.
Also, GRASS already includes allocation algorithms.

We could also imagine performing shortest path calculations for extra
vertexes that aren't initally on the graph: shop locations, particular
events, houses...

For this, spatstat::project2segment might be the appropriate function.

Based on your initial code, I think we could write a function that would
combine a linear network with extra vertexes into a single graph. The
vertex labels could discriminate the intersection nodes from the thematic
ones.


2012/8/31 Barry Rowlingson 

> On Fri, Aug 31, 2012 at 11:58 AM, Rowlingson, Barry
>  wrote:
>
> >  If I didn't have to be on a train in three hours I'd code this up...
>
>  Oh who am I kidding:
>
> buildTopo <- function(lines){
>   require(rgeos)
>   require(igraph)
>
>   g = gIntersection(lines,lines)
>   edges = do.call(rbind,lapply(g@lines
> [[1]]@Lines,function(ls){as.vector(t(ls@coords))}))
>   lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
>
>   froms = paste(edges[,1],edges[,2])
>   tos = paste(edges[,3],edges[,4])
>
>   graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
>   E(graph)$weight=lengths
>   return(graph)
>
> }
>
>  > m = readOGR(".","Line")
>  > gg = buildTopo(m)
>
> compute shortest distance between two vertices (note the -1 is because
> I'm still on old igraph with 0-indexing):
>
>  > E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]
>
> Edge sequence:
>
> [0]  0.81711193 47.57711634 -- -0.28546154 48.4287593
> [1]  2.2151139 46.49728033  -- 0.81711193 47.57711634
> [9]  1.57160852 45.36348513 -- 2.2151139 46.49728033
> [10] 1.39973834 45.06066625 -- 1.57160852 45.36348513
> [11] 1.29755246 44.88062446 -- 1.39973834 45.06066625
> [12] 0.91544563 43.67971729 -- 1.29755246 44.88062446
> [13] 0.88815229 43.43407719 -- 0.91544563 43.67971729
>
> I think this could be improved by adding the coordinates to the nodes.
> But anyway, 90% done.
>
> Barry
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Lines Network to Graph

[Mathieu Rajerison]

Hi,


Considering a linear network, how can I generate a graph of it that I could
use in igraph? So that I could perform shortest paths and other
calculations...

The vertexes would be the intersections of lines.


Thanks

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Create raster based on overlapping

[Mathieu Rajerison]

See raster::mask and raster::cover

2012/8/8 Thiago Veloso 

>  Dear all,
>
>  I have two rasters: r1 is a 0.5 degree netcdf file produced by a
> vegetation model. r2 is a 30m shapefile which I rasterized to a raster with
> the same 30m resolution. r2 values are all NA and 1.
>
>  Please see below:
>
> > r1
> class   : RasterLayer
> dimensions  : 10, 18, 180  (nrow, ncol, ncell)
> resolution  : 0.5, 0.5  (x, y)
> extent  : -53, -44, -25, -20  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> values  : in memory
> layer name  : variable
> zvar: variable
>
> > r2
> class   : RasterLayer
> dimensions  : 14469, 22190, 321067110  (nrow, ncol, ncell)
> resolution  : 0.0002702699, 0.0002702702  (x, y)
> extent  : -52.60649, -46.6092, -23.7172, -19.80666  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
> values  : in memory
> layer name  : cs_2005spx1
> min value   : 1
> max value   : 1
>
>  Here you can see a join plot:
> http://img855.imageshack.us/img855/2820/plotzoom.png
>
>  Is it possible to create a new raster, which is composed only by r1 cells
> which are "touched" by r2 cells?
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Google Mercator EPSG:900913, in which init file to add it?

[Mathieu Rajerison]

Hi,

I'm using Google Mercator (EPSG:900913) for one purpose but it doesn't seem
recognized by rgdal..

Where to add it please?

Thanks!

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Distance

[Mathieu Rajerison]

Hi,


>From scratch and using spatstat package, I'd say yes.

Convert your polygons into psp objects

Using angles, generate for each point a line passing through it and with
the extent of your data.

Find crossing points and calculate the distance between your point and each
corresponding crossing point. Find the minimum distance and the point.

Some functions: psp and crossing.psp


Best,

Mathieu

2012/7/6 Jocelyn Esquivel 

> Hi useRs
>
> I need to calculate the near distance between points (shapefile 1) and a
> polygon border (shapefile 2). I need to set a specific angle to calculate
> the linear distance from each point to the border of the polygon, but
> considering a specific direction. Is this possible?
>
> Thanks!!
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Conflict of raster package with function SMA of package TTR

[Mathieu Rajerison]

Hi,

To specify from which package the function has to be used

TTR::SMA

2012/7/6 Klaus Jacobi 

> Robert,
>
> as a hydrologist I am often calculating moving averages. For that,
> function SMA of package TTR was very convenient when dealing with time
> series.
>
> Yesterday I noted, that applying the function gave a strange error:
>
> > Ts.df[,4] <- SMA(Ts.df[,3],n=i) # moving average, at end of period
> Error in function (classes, fdef, mtable)  :
>   unable to find an inherited method for function "reclass", for signature
> "numeric"
>
> If some time of fiddling around with the data, I noticed, that the error
> occurs only when package "raster" is loaded. As soon as I detach it, the
> command is executed correctly.
>
> When I load library "TTR" first and afterwards "raster", I get the
> following messages:
>
> Loading required package: xts
> Loading required package: zoo
>
> Attaching package: ‘zoo’
>
> The following object(s) are masked from ‘package:base’:
>
> as.Date, as.Date.numeric
>
>
> Attaching package: ‘xts’
>
> The following object(s) are masked from ‘package:raster’:
>
> reclass
>
> I can even have raster and TTR simultaneously loaded, as long as I load
> TTR first and than raster. The other sequence (I am using package "raster"
> more frequently than TTR) causes the problem.
>
> Is there a change to fix that or should I just live with this workaround
> (i.e. loading packages in the right sequence).
>
> Klaus Jacobi
> k_jac...@surf2000.de
> www.disk-world.net
>
>
>
>
> [[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Species richness in R

[Mathieu Rajerison]

Hi,


You could use over( ) function in order to compute the number of different
species.

You could also use spatstat::marktable

Some ideas here:
https://stat.ethz.ch/pipermail/r-sig-geo/2011-May/011722.html

In vegan package, you can find a Simpson and Shannon diversity calculation
function: vegan::diversity


2012/7/2 Barnabas Daru 

> Dear All,
>
> I write to ask if anyone has an R code to calculate species richness per
> quarter degree grid.
>
>
>
> I have a set of polygon shapefiles showing distribution of different
> species within a country.
>
> I would like to first create quarter degree grids over the entire country
> and using the polygon shapefiles, to extract species richness per grid.
>
>
>
> Any help with R code would be highly appreciated.
>
> Barnabas Daru
>
>
>   \-/
>/\
>   /--|
>  /---/ Barnabas Daru
>  |--/  PhD Student,
>  \-/   African Centre for DNA Barcoding,
>  /\University of Johannesburg,
> /--\   PO Box 524, Auckland Park, 2006,
> |---\  South Africa.
>  \---\ Lab: +27 11 559 3477
>   \--| Mobile: +277 3818 9583
>\-/ Website: http://acdb.co.za/index.php?page=mr-barnabas-h-daru
>/\
>   /--\
>
> #…if you can think it, you can do it.
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] breaking lines at intersection but without GRASS

[Mathieu Rajerison]

Hi,


My final goal is to generate a graph from a road network. I didn't manage
to install shp2graph correctly, so I want to do it from scratch.

Is there a way to break lines at intersection but without rgrass?

Also, apart from shp2graph which seems a very good package, does one exist
that would turn my road network into a graph?


Thanks!

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] External "Envelope" of a set of points in space - not convex hull

[Mathieu Rajerison]

Hi,

alphahull is nice but having a SpatialPolygonsDataFrame from the result is
quite tricky

If anyone has the solution..

2012/6/22 Michael Sumner 

> Try the alphahull package, which is on CRAN.
>
> I don't believe it will smooth the result for you, but you could take
> it from there and use other R functions to do that.
>
> Cheers, Mike.
>
>
> On Fri, Jun 22, 2012 at 4:56 PM, PCB  wrote:
> > Dear All,
> >
> > I have a set of observations on an area, and I would like to define the
> > "investigated area" as the smallest polygon (preferably smoothed with
> e.g.
> > splines fit) that encloses all the points. This is similar to the convex
> > hull problem, except that my data have "gulfs", i.e. areas that are not
> > sampled (and should not be), between two "fingers" of areas
> investigated. I
> > am sure I have seen something about solving this lately, but after a
> couple
> > days searching could not find it (I am sure I am using the wrong
> > keywords...).
> >
> > Any help would be appreciated
> >
> > --
> > View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/External-Envelope-of-a-set-of-points-in-space-not-convex-hull-tp7580290.html
> > Sent from the R-sig-geo mailing list archive at Nabble.com.
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> Michael Sumner
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Help generating SpatialGridDataFrame from matrix

[Mathieu Rajerison]

Hi,

Maybe you could use raster::resample

2012/6/18 Jordan Winkler 

> I am confused on how to generate spatial grid data frames with attribute
> data and I've come across many different methods.  I am ultimately trying
> to extract temperature data from the CRU 3.1 data set for cropped areas
> globally.  I have both the CRU data and the cropped area dataset (Monfreda)
> in NetCDF.  I have the lon/lat vectors as well as a matrix of 1 timestep of
> the temperature data, and the same for the cropped area.  The main issue is
> that these two sets are of different resolution.  The crop data are 5
> minute, the temp is 1 degree.  My thought was to either extract the
> temperature values at each cell in the crop data set with a value > 0, then
> extract this to the country level using the appropriate shapefile.  Any
> guidance on where to start would be much appreciated.
>
> --
> --
> Jordan R Winkler
> Teaching Fellow / PhD Candidate
> Dept. of Geography and Environment
> Boston University
> 675 Commonwealth Avenue, STO 338
> Boston, MA 02215
>
>[[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] spatiallines / psp object to igraph networks

[Mathieu Rajerison]

Hi,


I've found the shp2grpah package which transforms a Lines network to an
object that can be used within igraph.

But do you know another package / other methods to coerce spatiallines/psp
objects into graphs?


Thanks

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Making a raster from a shapefile using function rasterize with a categoric field

[Mathieu Rajerison]

Hi,

It surely is...

You could transform lu$COB_AGROP as numeric like this:
lu$COB_AGROP_num <- as.numeric(as.factor(lu$COB_AGROP))

lur =  rasterize(lu, crr, field = "COB_AGROP_num")

Best,

Mathieu

2012/5/24 Manuel Spínola 

> Dear list members,
>
> I am trying to make a raster from a shapefile that I imported to r using
> the function "readOGR".  When I use the function rasterize and try to plot
> the object nothing is plotted.   Is this happening because the field is
> categoric (factor) and not numeric?
>
> > cr =
>
> readOGR(dsn="/Users/manuelspinola/Documents/ProyectosRespacial/EvalHab_reporte",
> layer="Cr_wgs84_meso")
>
> > crproy = spTransform(cr, CRS("+proj=tmerc +lon_0=-84 +lat_0=0 +x_0=50
> +k=0. +datum=WGS84"))
>
> > lu =
>
> readOGR(dsn="/Users/manuelspinola/Documents/ProyectosRespacial/Ecologia_del_paisaje",
> layer="usosuelos2006_cr05")
>
> > crr = raster(crproy)
>
> > lur =  rasterize(lu, crr, field = "COB_AGROP")
> Found 157480 region(s) and 227521 polygon(s)
> Mensajes de aviso perdidos
> 1: In .polygonsToRaster(x, y, ...) : selected field is character type
> 2: In .polygonsToRaster(x, y, ...) : NAs introducidos por coerción
>
> > lur
> class   : RasterLayer
> dimensions  : 10, 10, 100  (nrow, ncol, ncell)
> resolution  : 37190.08, 35181.17  (x, y)
> extent  : 286586.4, 658487.2, 889526.4, 1241338  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=tmerc +lon_0=-84 +lat_0=0 +x_0=50 +k=0.
> +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> values  : in memory
> min value   : NA
> max value   : NA
>
> plot(lur)
>
> Best,
>
> Manuel
>
>
> --
> *Manuel Spínola, Ph.D.*
> Instituto Internacional en Conservación y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspin...@una.ac.cr
> mspinol...@gmail.com
> Teléfono: (506) 2277-3598
> Fax: (506) 2237-7036
> Personal website: Lobito de río <
> https://sites.google.com/site/lobitoderio/>
> Institutional website: ICOMVIS 
>
>[[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] plot.raster - rgb(0,0,0,x) - jenks intervals

[Mathieu Rajerison]

Hi,

Finally, I found how to do it.

Using break argument in plot function.

2012/5/23 Mathieu Rajerison 

> Hi,
>
>
> I have a raster r that 'd like to superimpose on another raster with
> adjusted alpha values.
>
> The color would be black
>
> If I do:
> pal <- sapply(0:100, function(x) rgb(0,0,0,x/100))
> plot(r, col=pal)
>
> , my palette is adjusted with alpha values by equal interval distribution.
>
> Is there a way to plot the raster with another classification, by example
> jenks?
>
>
> Any help would be greatly appreciated,
>
> Mathieu
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] plot.raster - rgb(0,0,0,x) - jenks intervals

[Mathieu Rajerison]

Hi,


I have a raster r that 'd like to superimpose on another raster with
adjusted alpha values.

The color would be black

If I do:
pal <- sapply(0:100, function(x) rgb(0,0,0,x/100))
plot(r, col=pal)

, my palette is adjusted with alpha values by equal interval distribution.

Is there a way to plot the raster with another classification, by example
jenks?


Any help would be greatly appreciated,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] image segmentation in R - forest stand delineation from LiDAR data

[Mathieu Rajerison]

Hi!

Maybe some ideas here:
http://casoilresource.lawr.ucdavis.edu/drupal/node/993 using randomForest ?

2012/5/20 Omphalodes Verna 

> Dear list members!
>
> I kindly ask you for advice in image segmentation in R. I have read a few
> posts in R-sig-Geo list, but I did not found clear answer. I tried k-means
> clustering (in package ''raster''), but my input images is very
> heterogeneous (forest stand parameters, derived form LiDAR data in uneven
> aged, mixed forest, with small scale disturbances) and output from
> clustering process is very fine scale image with high local variation. I am
> looking for image segmentation for forest stand delineation.
>
> Thanks to all for answers.
>
> OV
>
>[[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] plotting raster - color values per pixel

[Mathieu Rajerison]

Hi!

Finally i chose the plotRGB method.

Here is the code:
cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
cols <- col2rgb(cols)

red <- r; values(red) <- cols[1, ]
green <- r; values(green) <- cols[2, ]
blue <- r; values(blue) <- cols[3, ]

r.rgb <- stack(red, green, blue)

plotRGB(r.rgb)

Thanks!

2012/5/11 Barry Rowlingson 

> On Fri, May 11, 2012 at 9:20 AM, Mathieu Rajerison
>  wrote:
>
> >> In col, I must specify a palette, but how can I use my color vector?
> >>
> >> Here is the line of code for the moment (which doesn't work):
> >> cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
> >> plot(r, col=topo.colors(20), axes=FALSE)
>
>  You can do this by instead of plotting the 'r' raster (which is your
> values) you plot a raster containing the values 1 to nrows*ncols. and
> pass col as a vector of colours of length nrows*ncols. Then each cell
> in the raster maps exactly to its colour value.
>
>  Disadvantage: the legend is now 1:(nrows*ncols)
>
>  The other way is to create a 3-band raster from the red, green, blue
> components of the colours.
>
> Barry
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] plotting raster - color values per pixel

[Mathieu Rajerison]

Excuse me, I've made a little mistake:
cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
plot(r, col=*cols*, axes=FALSE)

2012/5/11 Mathieu Rajerison 

> Hi,
>
>
> I'm trying to get an isarithmic map of votes in France.
>
> I've calculated a vector of n colors depending on the pixel value, n being
> the raster number of pixels.
>
> In col, I must specify a palette, but how can I use my color vector?
>
> Here is the line of code for the moment (which doesn't work):
> cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
> plot(r, col=topo.colors(20), axes=FALSE)
>
>
> Thanks!
>
> Mathieu
>
> --
> the entire code:
> library(raster)
> library(gstat)
> x <- g2.pp$x; y <- g2.pp$y
> z <- marks(g2.pp)
> df <- as.data.frame(cbind(x, y, z))
> r <- raster(extent(g2))
> res(r) <- 500
> mg.pchol <- gstat(id="z", formula=z~1, locations=~x+y, data=df, nmax=9,
> set=list(idp=.5))
> r <- interpolate(r, mg)
> cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
> plot(r, col=topo.colors(20), axes=FALSE)
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] plotting raster - color values per pixel

[Mathieu Rajerison]

Hi,


I'm trying to get an isarithmic map of votes in France.

I've calculated a vector of n colors depending on the pixel value, n being
the raster number of pixels.

In col, I must specify a palette, but how can I use my color vector?

Here is the line of code for the moment (which doesn't work):
cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
plot(r, col=topo.colors(20), axes=FALSE)


Thanks!

Mathieu

--
the entire code:
library(raster)
library(gstat)
x <- g2.pp$x; y <- g2.pp$y
z <- marks(g2.pp)
df <- as.data.frame(cbind(x, y, z))
r <- raster(extent(g2))
res(r) <- 500
mg.pchol <- gstat(id="z", formula=z~1, locations=~x+y, data=df, nmax=9,
set=list(idp=.5))
r <- interpolate(r, mg)
cols <- hsv(ifelse(values(r) > .5, 0, .7), values(r), 1)
plot(r, col=topo.colors(20), axes=FALSE)

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] analyzing spatial segregation of sexes

[Mathieu Rajerison]

Hi!

Why not use, for example, KCross spatstat function, given two point
patterns: one for male individuals and one for female ones?


2012/5/9 Wall, Wade A ERDC-RDE-CERL-IL 

> Hi all,
>
> I have a data set that consists of x,y coordinates for individual plants,
> along with sex of the individual.
>
> Example
>
> Ind sex x   y
> 1   F   10  5
> 2   M   9   4
> . . .
> 200 M   20  4
>
> for 20 populations of the species. I would like to know if males and
> females are randomly dispersed, or tend to segregate.
>
> I would like to use Philip Dixon's method (1994), and there is an
> implementation in the R package "spatialsegregation" (dixon(X, prepR=0)).
> However,
> I am not very familiar with spatial analyses in R and am not sure how to
> structure the input object X. The function description says "Bivariate i.e.
> 2-type point pattern (see package 'spatstat')"
>
> Does anyone know how to structure the data that I have (example above) so
> that I can use dixon()? Thanks for any help.
>
> Wade A. Wall
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Visibility Analysis in R

[Mathieu Rajerison]

Hi,

I had written a post about intervisibility with R & GRASS on my blog. I
hope it will be helpful for you. I don't explain how to connect R and GRASS
but you will find a lot of articles about that on the web.
http://www.datagistips.com/2011/09/on-road-with-r-grass-intervisibility.html

Best

2012/4/21 Barry Rowlingson 

> On Sat, Apr 21, 2012 at 1:38 PM, Raphael Fuhrer
>  wrote:
> > Dear all
> >
> > Does anybody know if it is possible to conduct a visibility analysis in
> R?
> > I have
> > - a large set of points = observer points
> > - a digital surface model as an ascii grid
> > - a set of points and polygons defining an amount of grid cells of
> interest
> > For every individual observer point, the task is to check how many cells
> of interest are visible.
> >
> > Since the number of observer points is limited to 8 in ArcGIS and
> calculation time is far too long, I decided to use R. Unfortunately, I
> could not find any helpful packages or information on the internet.
>
>  I suspect your best option is to use GRASS-GIS, and there is a
> package for R-GRASS integration but I don't know much about it:
>
> http://cran.r-project.org/web/packages/spgrass6/index.html
>
> The relevant GRASS module is r.viewshed (r here stands for raster, not R):
>
> http://grass.osgeo.org/manuals/html70_user/r.viewshed.html
>
> Free and open source, of course it is!
>
> Barry
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spdep::listw2lines - removing duplicated lines

[Mathieu Rajerison]

Hello Roger,

Thank you a lot for the answer.

The goal was just creating a map using R and Processing to represent
relationships as a neuronal network. (just having fun :))
http://www.flickr.com/photos/10519370@N04/6921277544/

Mathieu


2012/4/11 Roger Bivand 

> On Wed, 11 Apr 2012, Mathieu Rajerison wrote:
>
>  Hi,
>>
>>
>> I'm using spdep to generate neighbors, than listw2lines to create a
>> SpatialLinesDataFrame.
>>
>> Considering two neighbors N1 and N2, the line will be duplicated: one
>> N1->N2 and one N2->N1
>>
>
> Only if the neighbours are symmetric, and in graph-based or k-nearest
> neighbours, this is not the case.
>
>
>
>> Would you have any idea on how to select only one line from the two, for
>> example the first?
>>
>
> If the neighbours are symmetric, it would be possible to do for example:
>
> library(spdep)
> set.seed(1)
> xy <- cbind(x=runif(50), y=runif(50))
> nb <- dnearneigh(xy, 0, 0.2)
> nb
> is.symmetric.nb(nb)
> nb1 <- vector(mode="list", length=length(nb))
> for (i in seq(along=nb)) nb1[[i]] <- nb[[i]][nb[[i]] <= i]
> class(nb1) <- "nb"
> attr(nb1, "region.id") <- attr(nb, "region.id")
> nb1
>
> It could probably also be done with gEquals() in rgeos also for asymmetric
> neighbours to remove duplications where they occur.
>
> Roger
>
>
>  Actually, it would consist in removing duplicated lines with same geometry
>> but different directions.
>>
>>
>> Thanks in advance, any help would be greatly appreciated!
>>
>> Mathieu
>>
>>[[alternative HTML version deleted]]
>>
>> __**_
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>
>>
> --
> Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: roger.biv...@nhh.no
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] spdep::listw2lines - removing duplicated lines

[Mathieu Rajerison]

Hi,


I'm using spdep to generate neighbors, than listw2lines to create a
SpatialLinesDataFrame.

Considering two neighbors N1 and N2, the line will be duplicated: one
N1->N2 and one N2->N1

Would you have any idea on how to select only one line from the two, for
example the first?
Actually, it would consist in removing duplicated lines with same geometry
but different directions.


Thanks in advance, any help would be greatly appreciated!

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] buffer between points and line shapefiles

[Mathieu Rajerison]

Hi,

Maybe
spatstat::nncross?


2012/4/9 Juan Tomas Sayago 

> Hello
> This question has been probably addressed by someone else and responded but
> could not find it in the archive so, here it is.
> I am trying to create a buffer of distance around points and then check if
> the point is located within some distance of a river (line shapefile). Does
> anyone has idea about it.
> Thanks
> Juan
>
> --
> Juan Tomás Sayago Gómez
> Graduate Research Assistant
> West Virginia University - RRI
> 886 Chestnut Ridge Road, Room 520
> P.O. Box 6825
> Morgantown, WV 26506-6825
>
>[[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] One lot <-> several owners : how to duplicate geometries?

[Mathieu Rajerison]

Here is another method using rgeos..
The polygon geometry will be stored in a column and re-generated after the
data mapping is done

spP$WKT <- writeWKT(spP, byid=TRUE)

# I merge
spP.df <- as.data.frame(spP)
m <- merge(as.data.frame(spP), df, by.x="IDLOT", by.y="IDLOT")

# row names for m
IDs <- as.character(m$IDLOT)
for (i in levels(factor(IDs)))
{
  loc <- which(IDs==i); len = length(loc)
  IDs[loc] <- sapply(1:len, function(j) paste(IDs[loc][j], j, sep="."))
}
rownames(m) <- IDs

# SpatialPolygonsDataFrame
out <- vector(mode="list")
for (i in 1:nrow(m)) {
  out[[i]] <- readWKT(m$WKT[i], id = IDs[i], p4s = CRS("+init=EPSG:2154"))
}
save(out, file="out.RData")
r <- do.call("rbind", out)
spPdf <- SpatialPolygonsDataFrame(r, data = m)


All this is quite complex to do in R, but easier with PostgreSQL, except if
another method, simpler, exists?

Best,

Mathieu


2012/3/14 Mathieu Rajerison 

> Hi,
>
> Finally, I managed to do it. Here is how I processed. Maybe there was a
> simpler way. I had to build a SpatialPolygons from scratch..
>
> # On one side, I had a SpatialPolygonsDataFrame called spP
> # One the other side, a data frame called df
> spP <- readOGR(".", "LOTS")
>
> # I create a column IDGEOM:
> spP$IDGEOM <- 1:nrow(spP)
>
> # I merge spP and df
> spP.df <- as.data.frame(spP)
> m <- merge(spP.df,df, by.x="IDLOT", by.y="IDLOT")
>
> # Then, I give each line a unique id based on a numbered IDLOT
> # If the IDLOT 1320 appears twice (two owners for the lot), it will appear
> like this: 1320.1 then 1320.2
> # These IDs will be the Polygon IDs when I will build the SpatialPolygons
> object
> IDs <- as.character(m$IDLOT)
> for (i in levels(factor(IDs)))
> {
>   loc <- which(IDs==i); len = length(loc)
>   IDs[loc] <- sapply(1:len, function(j) paste(IDs[loc][j], j, sep="."))
> }
> rownames(m) <- IDs
>
> # Then, I create a Polygons List and finally the SpatialPolygonsDataFrame
> pl <- lapply(1:nrow(m), function(i) Polygons(slot(slot(spP,
> "polygons")[[m$IDGEOM[i]]], "Polygons"), ID=IDs[i]))
> spPdf <- SpatialPolygonsDataFrame(SpatialPolygons(pl), data = m)
>
>
> If anyone has a better idea or can make the code cleaner..
> There may be a way to do it a similar way with rgeos and WKT (?)
>
> Best,
>
> Mathieu
>
>
> 2012/3/14 Mathieu Rajerison 
>
>> Hi,
>>
>> I've already sent this e-mail to r-sig-...@stat.math.ethz.ch, but the
>> correct email address may be r-sig-geo@r-project.org
>>
>> I have on one side geometries of lots that are unique
>> (spatialPolygonsDataFrame). On the other side, I have a huge file of owners
>> affected to each lot (data.frame).
>>
>> I know how to join one geometry to data with merge function, keeping row
>> names equivalent between geodata and data but here, as one lot can have
>> many owners, I need to duplicate the geometries when I meet many owners for
>> one lot.
>>
>> Is there a way to do so?
>>
>> Thanks in advance,
>>
>> Mathieu
>>
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] One lot <-> several owners : how to duplicate geometries?

[Mathieu Rajerison]

Hi,

Finally, I managed to do it. Here is how I processed. Maybe there was a
simpler way. I had to build a SpatialPolygons from scratch..

# On one side, I had a SpatialPolygonsDataFrame called spP
# One the other side, a data frame called df
spP <- readOGR(".", "LOTS")

# I create a column IDGEOM:
spP$IDGEOM <- 1:nrow(spP)

# I merge spP and df
spP.df <- as.data.frame(spP)
m <- merge(spP.df,df, by.x="IDLOT", by.y="IDLOT")

# Then, I give each line a unique id based on a numbered IDLOT
# If the IDLOT 1320 appears twice (two owners for the lot), it will appear
like this: 1320.1 then 1320.2
# These IDs will be the Polygon IDs when I will build the SpatialPolygons
object
IDs <- as.character(m$IDLOT)
for (i in levels(factor(IDs)))
{
  loc <- which(IDs==i); len = length(loc)
  IDs[loc] <- sapply(1:len, function(j) paste(IDs[loc][j], j, sep="."))
}
rownames(m) <- IDs

# Then, I create a Polygons List and finally the SpatialPolygonsDataFrame
pl <- lapply(1:nrow(m), function(i) Polygons(slot(slot(spP,
"polygons")[[m$IDGEOM[i]]], "Polygons"), ID=IDs[i]))
spPdf <- SpatialPolygonsDataFrame(SpatialPolygons(pl), data = m)


If anyone has a better idea or can make the code cleaner..
There may be a way to do it a similar way with rgeos and WKT (?)

Best,

Mathieu

2012/3/14 Mathieu Rajerison 

> Hi,
>
> I've already sent this e-mail to r-sig-...@stat.math.ethz.ch, but the
> correct email address may be r-sig-geo@r-project.org
>
> I have on one side geometries of lots that are unique
> (spatialPolygonsDataFrame). On the other side, I have a huge file of owners
> affected to each lot (data.frame).
>
> I know how to join one geometry to data with merge function, keeping row
> names equivalent between geodata and data but here, as one lot can have
> many owners, I need to duplicate the geometries when I meet many owners for
> one lot.
>
> Is there a way to do so?
>
> Thanks in advance,
>
> Mathieu
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] One lot <-> several owners : how to duplicate geometries?

[Mathieu Rajerison]

Hi,

I've already sent this e-mail to r-sig-...@stat.math.ethz.ch, but the
correct email address may be r-sig-geo@r-project.org

I have on one side geometries of lots that are unique
(spatialPolygonsDataFrame). On the other side, I have a huge file of owners
affected to each lot (data.frame).

I know how to join one geometry to data with merge function, keeping row
names equivalent between geodata and data but here, as one lot can have
many owners, I need to duplicate the geometries when I meet many owners for
one lot.

Is there a way to do so?

Thanks in advance,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] One lot <-> several owners : duplicating geometries

[Mathieu Rajerison]

Hi,

I have on one side geometries of lots that are unique
(spatialPolygonsDataFrame). On the other side, I have a huge file of owners
affected to each lot (data.frame).

I know how to join one geometry to data with merge function, keeping row
names equivalent between geodata and data but here, as one lot can have
many owners, I need to duplicate the geometries when I meet many owners for
one lot.

Is there a way to do so?

Thanks in advance,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] tilted satellite projection

[Mathieu Rajerison]

Hi,

Does anybody know if we can use, in R, the tilted satellite projection?
http://user.gs.rmit.edu.au/rod/files/publications/Tilted%20Camera%20Map%20Projection.pdf

Best,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] dismo::gmap: can not create a RasterLayer from this file

[Mathieu Rajerison]

Robert,

Maybe it's because I'm behind a proxy?

2012/2/9 Robert J. Hijmans 

> Mathieu,
> This works for me. Apparently something went wrong with the download:
> "taille téléchargée 270471 != taille déclarée 200". But why?
> Robert
>
> On Tue, Feb 7, 2012 at 4:39 AM, Mathieu Rajerison <
> mathieu.rajeri...@gmail.com> wrote:
>
>> Hi,
>>
>>
>> I'm trying to use gmap function from dismo.
>>
>> I apply the example but it gives me an error message 'Cannot create a
>> RasterLayer object from this file.'
>>
>> > library(dismo)
>> > e = extent( -121.9531 , -120.3897 , 35.36 , 36.61956 )
>> > r = gmap(e)
>> Erreur dans .rasterObjectFromFile(x, band = band, objecttype =
>> "RasterLayer", :
>> Cannot create a RasterLayer object from this file.
>> De plus : Message d'avis :
>> In download.file(gurl, filename, mode = "wb", quiet = TRUE) :
>> taille téléchargée 270471 != taille déclarée 200
>>
>> Any help would be greatly appreciated!
>>
>>
>> Thanks,
>>
>> Mathieu
>>
>>[[alternative HTML version deleted]]
>>
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] dismo::gmap: can not create a RasterLayer from this file

[Mathieu Rajerison]

Hi,


I'm trying to use gmap function from dismo.

I apply the example but it gives me an error message 'Cannot create a
RasterLayer object from this file.'

> library(dismo)
> e = extent( -121.9531 , -120.3897 , 35.36 , 36.61956 )
> r = gmap(e)
Erreur dans .rasterObjectFromFile(x, band = band, objecttype =
"RasterLayer", :
Cannot create a RasterLayer object from this file.
De plus : Message d'avis :
In download.file(gurl, filename, mode = "wb", quiet = TRUE) :
taille téléchargée 270471 != taille déclarée 200

Any help would be greatly appreciated!


Thanks,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] finding coordinates of points on a dilation contour

[Mathieu Rajerison]

Hi,


Let's say your contour initially is an owin object named W

you apply
Wd <- dilation(W, 10)

Then, you coerce your window to a psp object Sg
Sg <- as.psp(Xd)

Finally, you can get the coordinates of the points with
Sg$ends

If your request consists in getting uniformally distributed points along
the dilated contour, you could use the functions pointsOnLines or
runifpointOnLines


Best,

Mathieu

2012/1/31 Lorenzo Cattarino 

> Hi,
>
> I have a set of coordinates representing a shape (i.e. a square). I would
> like to perform an Euclidean dilation operation, i.e., find all the point
> that lie at a distance r from the shape contour. I tried with the dilation
> function in the "spatstat" package but I could not get the coordinates of
> the points forming the dilated contour, nor the number of those points. In
> the case of a square it might be a trivial exercise. However I don't always
> have squared objects to perform the dilation on.
>
> coord <-
> matrix(c(1,1,1,2,1,3,1,4,1,5,1,6,2,1,3,1,4,1,5,1,6,1,2,6,3,6,4,6,5,6,6,6,6,2,6,3,6,4,6,5),ncol=2,byrow=T)
>
> I would really appreciate if anyone can point me in the right direction.
>
> Cheers
> Lorenzo
>
>
>[[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] R package to know distances between city [or zipcode]

[Mathieu Rajerison]

Hi,

You could use spatstat::nndist/nnwhich to find the nearest points.

Best,

Mathieu

2012/1/27 Antoine Lucas 

> Dear All,
>
> I have build a small package, that main object is to compute:
>
> 1. city names that are close to another city
> 2. city names that are close to a zipcode
>
> I have two questions:
>
> -
> a. for now, dataset is french data (zipcode + city + lat,lon): I would
> like someting like:
>
> library("mynewpackage")
>
> setDefaultCountry("fr")
> getCityNearToZipCode(1234)  => french city list
>
> setDefaultCountry("uk")
> getCityNearToZipCode(1234)  => uk city list
>
> for now, there is only one dataset that I use directly in my R functions...
>
> -
> b. Is this package needed, or is it redunctdant with an existing package ?
>
> Antoine.
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] cut a line at specific point

[Mathieu Rajerison]

Hi,

Do you have one layer of segments or two?

You'll find what you need in the spatstat package.

See
(self)crossing.psp find crossing points between two line segment patterns
Afterwards, it must possible to build new segments based on the
intersection points

Also, you have the GRASS v.clean that allows breaking lines at intersection
points

Best,

Mathieu

2012/1/26 Michiel Faber 

> Dear all,
>
> I'am using R extensively privately and during work.
> I have just discovered the Spatial classes of R and want to use them in
> my home project. I want to find the intersections between lines and cut
> the lines at those intersections.
> For example: Two lines crossing each other once will result in four new
> lines.
>
> I have searched intensively, but I can't find a simple way of cutting a
> (Spatial) line in two. I can make a data.frame out of it, cut it in two
> and make two (Spatial) lines again, but that doesn't seem logical.
>
> Is there a simple way to do it?
> The gIntersection() between the lines in the example yields two
> intersections (int.s1_s2). I would like to have six new lines, each line
> being cut in three parts by the points in int.s1_s2.
>
> Thanks for any pointers,
> Michiel
>
> library(rgeos)
> x1 <- c(3,  4,  6,  8, 10)
> y1 <- c(0.23, 0.48, 0.45, 0.08, 0.18)
> x2 <- c(1,  2,  6,  7, 10)
> y2 <- c(0.27, 0.36, 0.16, 0.30, 0.85)
>
> #is there also a easier way of doing this?
> s1 <- SpatialLines(list(Lines(list(Line(cbind(x1,y1))),ID=1)))
> s2 <- SpatialLines(list(Lines(list(Line(cbind(x2,y2))),ID=1)))
>
> int.s1_s2 <- gIntersection(s1,s2)
>
> plot(s1, type="b", col="green", xlim=c(1,10))
> lines(s2, type="b", col="blue")
> points(int.s1_s2, pch=4, col="black")
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Lenght of Common border

[Mathieu Rajerison]

Hi,

This function is available in spgrass::vect2neigh

"vect2neigh return area neighbours with shared boundary length"

It requires installing GRASS


2012/1/23 Gaston Pezzuchi 

> Hi there:
>
> I am trying to find the actual formula to calculate the length of the
> common border between adjacent  polygons. So far I've been using:
>
> - Let r1 and r2 be the perimeters of Polygons 1 and 2.
> - Make R = r1 + r2
> - Make r = perimeter of the union between r1 and r2
> - Length of common border = (R - r)/2
>
>
>   This seems to work well if the common border is a single line segment,
> yet
> if the common border is comprised by more than one segment (two or more
> disjoint line segments), it will not work.
>
>Any help would be greatly appreciated.
>
> Regards
>
> Gaston
>
> Lic. Gastón Pezzuchi
> Magister UNIGIS en Sistemas de Información Geográfica
>
>[[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] extract value of polygone

[Mathieu Rajerison]

Hi,

First, you should try to get a continuous topography raster map from your
contour map.

Then, you could extract values at some points using raster::extract.

Best,

Mathieu

2012/1/18 Saman Monfared 

> Hi.
> I try to extract values of temperature in a topography  map. It is
> like a contour plot.
> I need the coordinates and temperature of some points.
> how can I perform it??
> thanks.
>
>
> --
>
> Saman Monfared
> Msc Student,
> Department of Statistics, Shiraz University,
> Shiraz 71454, Iran
> Email: samanmonfar...@gmail.com
>
> Tel: +98 917 5305167
>
>
>
> --
>
> Saman Monfared
> Msc Student,
> Department of Statistics, Shiraz University,
> Shiraz 71454, Iran
> Email: samanmonfar...@gmail.com
>
> Tel: +98 917 5305167
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] NYT map with non-overlapping labels

[Mathieu Rajerison]

Hi,


I've juste read this article about a map published on NYT
http://blog.revolutionanalytics.com/2012/01/nyt-uses-r-to-map-the-1.html

The map has been created using R. Actually, the map is all about
non-overlapping labels which algorithmic placement.
http://revolution-computing.typepad.com/.a/6a010534b1db25970b0168e5b1af04970c-pi

The article doesn't mention which function was used to perform this. It
just says maptools has been used

I just wondered if maptools::pointLabel function allows this kind of
process (?)


Best,

Mathieu

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo