Re: [R-sig-Geo] Spatial random forest prediction: Error when predicting at unseen locations at finer spatial scale
Hi Nikolaos:
As you can read in the github page of the package
(https://blasbenito.github.io/spatialRF/),
*"there are many things this package cannot do*:
*
Predict model results over raster data.
*
Predict a model result over another region with a different spatial
structure.
* ..."
More specifically, the error you report is probably related to a lot of
spatial predictors derived from the distance matrix (that are included
in the fitted model) that you are not passing to "predict" the model.
Cheers,
Marcelino
El 22/02/2023 a las 21:36, Nikolaos Tziokas escribió:
I am using the package spatialRF in R for a spatial random forest
regression (SRFR) task. I have one response variable and 4 predictors and I
am performing SRFR at a coarse spatial scale. My goal is to take the model
parameters and apply them to a finer spatial resolution in order to predict
the response variable at the finer spatial scale.
When I run
p <- stats::predict(object = model.spatial, #name of the
spatialRF model
data = s, # data.frame
containing the predictors at the fine spatial scale (without NaN values)
type = "response")$predictions
I am getting this error: Error in predict.ranger.forest(forest, data,
predict.all, num.trees, type,: Error: One or more independent variables not
found in data.
I have checked the column names of s and my original data.frame (the one I
used to build the model at the coarse scale) and they are the same. How can
I use the model i created at the coarse scale to predict the response
variable at a finer spatial scale?
Here is the code:
library(spatialRF)
library(stats)
wd = "path/"
block.data = read.csv(paste0(wd, "block.data.csv")) # coarse resolution
#names of the response variable and the predictors
dependent.variable.name <- "ntl"
predictor.variable.names <- colnames(block.data)[4:7]
#coordinates of the cases
xy <- block.data[, c("x", "y")]
block.data$x <- NULL
block.data$y <- NULL
#distance matrix
distance.matrix <- as.matrix(dist(block.data))
min(distance.matrix)
max(distance.matrix)
#distance thresholds (same units as distance_matrix)
distance.thresholds <- c(0, 20, 50, 100, 200, 500)
#random seed for reproducibility
random.seed <- 456
#creating and registering the cluster
local.cluster <- parallel::makeCluster(
parallel::detectCores() - 1,
type = "PSOCK")
doParallel::registerDoParallel(cl = local.cluster)
# fitting a non-spatial Random Forest
model.non.spatial <- spatialRF::rf(
data = block.data,
dependent.variable.name = dependent.variable.name,
predictor.variable.names = predictor.variable.names,
distance.matrix = distance.matrix,
distance.thresholds = distance.thresholds,
xy = xy,
seed = random.seed,
verbose = FALSE)
# Fitting a spatial model with rf_spatial()
model.spatial <- spatialRF::rf_spatial(
model = model.non.spatial,
method = "mem.moran.sequential",
verbose = FALSE,
seed = random.seed)
#stopping the cluster
parallel::stopCluster(cl = local.cluster)
# prediction at a finer spatial scale
s = read.csv(paste0(wd, "s.csv")) # df containg the predictors at fine
scale
p <- stats::predict(object = model.spatial,
data = s,
type = "response")$predictions
I tried solutions like:
levels(s$lc) <- levels(block.data$lc)
in case I had missing land cover types in the lc column between the spatial
scales, but it didn't work.
>From here
<https://drive.google.com/drive/folders/1KhnQEajpSKh59XuWkxTZcc_2YxPcxYW7?usp=sharing>
you can download the two data.frames.
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] st_intersection produce Geometry type: GEOMETRY instead of Geometry type: POLYGON
g_25000_c_polygons_only <- g_25000_c[st_is(x = g_25000_c, type = "POLYGON"),] #
with a comma
El 10/02/2023 a las 21:33, Manuel Spínola escribió:
Thnak you very much Bede-Fazekas.
I got the following error:
g_25000_c_polygons_only <- g_25000_c[st_is(x = g_25000_c, type = "POLYGON")]
Error in `[.data.frame`(x, i) : undefined columns selected
El vie, 10 feb 2023 a las 11:53, Bede-Fazekas Ákos ()
escribió:
Dear Manuel,
technically, the result of st_intersection(x, y), where both x and y are
POLYGONs, can be POINT, LINESTRING, POLGYON and GEOMETRY as well. The
result is GEOMETRY if the type of the different features is not the same
(e.g. POLYGON+POINT).
You can subset the result in this way:
g_25000_c_polygons_only <- g_25000_c[st_is(x = g_25000_c, type =
"POLYGON")]
HTH,
Ákos
___
Ákos Bede-Fazekas
Centre for Ecological Research, Hungary
2023.02.10. 18:32 keltezéssel, Manuel Spínola írta:
Dear list members,
I am trying to "crop" a polygon (grid) with a polygon, but the result is
an
sf object Geometry type: GEOMETRY, instead of an sf object Geometry type:
POLYGON.
How can I obtain an sf POLYGON?
nc = st_read(system.file("shape/nc.shp", package="sf"))
nc <- st_transform(nc, 3857)
g_5 <- st_make_grid(nc, cellsize = 5) |> st_as_sf()
g_5 <- g_5[nc, ]
g_5_d <- st_union(g_5)
g_25000 = st_make_grid(g_5_d, cellsize = 25000) |> st_as_sf()
g_25000 # Geometry type: POLYGON
g_25000_c <- st_intersection(g_25000, g_5_d)
g_25000_c # Geometry type: GEOMETRY
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Marcelino de la Cruz Rot
Coordinador funcional de Biología
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
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 nested grid in R
El 10/02/2023 a las 15:13, Roger Bivand escribió: On Fri, 10 Feb 2023, Marcelino de la Cruz Rot wrote: It depends on what you mean by "overimpose". Maybe this way? # some spatial polygon Please, in the spirit of the evolving r-spatial package ecosystem (maptools will retire during 2023): library(sf) R <- st_as_sf(affine(letterR, mat=diag(c(20,20 plot(R, border="blue") See https://r-spatial.org/r/2022/12/14/evolution2.html, maybe https://rsbivand.github.io/csds_jan23/bivand_csds_ssg_230117.pdf and https://www.youtube.com/watch?v=TlpjIqTPMCA=PLzREt6r1NenmWEidssmLm-VO_YmAh4pq9=1 Roger Thank you Roger! In the spirit of the evolving, I then suggest this to Manuel: # define this function (if it is not in your version of sf) st_as_sf.tess <- sf:::st_as_sf.owin # transform tesselations in sf's Dsf <- st_as_sf(D) Csf <- st_as_sf(C) Bsf <- st_as_sf(B) Asf <- st_as_sf(A) # crop tesselations with spatial polygon RDsf <- Dsf[R,] RCsf <- Csf[R,] RBsf <- Bsf[R,] RAsf <- Asf[R,] #overimpose cropped grids plot(RAsf) plot(RDsf, border="green", add=T) plot(RCsf, border="blue", add=T) plot(RBsf, border="red", add=T) plot(RAsf, add=T) Cheers, Marcelino library(maptools) R <- as(affine(letterR, mat=diag(c(20,20))), "SpatialPolygons") plot(R, border="blue") # "overimpose" grids A to C on R: plot(C, border="green", add=T) plot(B, border="red", add=T) plot(A, add=T) Cheers, Marcelino El 10/02/2023 a las 13:29, Manuel Spínola escribió: Thank you very much Marcelino. And how can overimpose those grids to a spatial polygon? Manuel On Fri, 10 Feb 2023 at 03:09 Marcelino de la Cruz Rot wrote: Dear Manuel, This is R. There is no "it is possible". Only "how" ;-). For example, with spatstat.geom, A <- tess(xgrid=seq(0,80, by=8),ygrid=seq(0,80, by=8)) B <- tess(xgrid=seq(0,80, by=4),ygrid=seq(0,80, by=4)) C <- tess(xgrid=seq(0,80, by=2),ygrid=seq(0,80, by=2)) D <- tess(xgrid=seq(0,80, by=1),ygrid=seq(0,80, by=1)) Cheers, Marcelino El 10/02/2023 a las 0:21, Manuel Spínola escribió: > Dear list members, > > Is it possible to generate a spatial nested grid in R? > > For example, a grid of several 8km x 8km tiles, and within that grid, I > want 4 tiles of 4km x 4km, and in each of those I want 4 tiles of 2km x > 2km, and in each of those I want 4 tiles of 1km x 1km. > > Manuel > -- Marcelino de la Cruz Rot Coordinador funcional de Biología Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo=05%7C01%7CRoger.Bivand%40nhh.no%7C0b29de36aca54242ddf008db0b6ade31%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638116325854200214%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Ngv4yEcnP6JM4zq5rpVgzcYgRHAFwPPS8JRFhMaRA58%3D=0 -- *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.cr <mailto:mspin...@una.ac.cr> mspinol...@gmail.com Teléfono: (506) 8706 - 4662 Institutional website: ICOMVIS <https://eur02.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.icomvis.una.ac.cr%2Findex.php%2Fmanuel=05%7C01%7CRoger.Bivand%40nhh.no%7C0b29de36aca54242ddf008db0b6ade31%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638116325854200214%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=fCKSKsYuj7JXD5CbsQDKFUkeAVpqekVl%2BbTYM2dpuzU%3D=0> Blog sobre Ciencia de Datos: https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmspinola-ciencia-de-datos.netlify.app%2F=05%7C01%7CRoger.Bivand%40nhh.no%7C0b29de36aca54242ddf008db0b6ade31%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638116325854200214%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=cVcjfP2hgnbfK%2BZCrCGFi4HhnQjGEWtx8zeJ%2BjBZDtE%3D=0 -- Marcelino de la Cruz Rot Coordinador funcional de Biología Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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 nested grid in R
It depends on what you mean by "overimpose".
Maybe this way?
# some spatial polygon
library(maptools)
R <- as(affine(letterR, mat=diag(c(20,20))), "SpatialPolygons")
plot(R, border="blue")
# "overimpose" grids A to C on R:
plot(C, border="green", add=T)
plot(B, border="red", add=T)
plot(A, add=T)
Cheers,
Marcelino
El 10/02/2023 a las 13:29, Manuel Spínola escribió:
Thank you very much Marcelino.
And how can overimpose those grids to a spatial polygon?
Manuel
On Fri, 10 Feb 2023 at 03:09 Marcelino de la Cruz Rot
wrote:
Dear Manuel,
This is R. There is no "it is possible". Only "how" ;-).
For example, with spatstat.geom,
A <- tess(xgrid=seq(0,80, by=8),ygrid=seq(0,80, by=8))
B <- tess(xgrid=seq(0,80, by=4),ygrid=seq(0,80, by=4))
C <- tess(xgrid=seq(0,80, by=2),ygrid=seq(0,80, by=2))
D <- tess(xgrid=seq(0,80, by=1),ygrid=seq(0,80, by=1))
Cheers,
Marcelino
El 10/02/2023 a las 0:21, Manuel Spínola escribió:
> Dear list members,
>
> Is it possible to generate a spatial nested grid in R?
>
> For example, a grid of several 8km x 8km tiles, and within that
grid, I
> want 4 tiles of 4km x 4km, and in each of those I want 4 tiles
of 2km x
> 2km, and in each of those I want 4 tiles of 1km x 1km.
>
> Manuel
>
--
Marcelino de la Cruz Rot
Coordinador funcional de Biología
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
*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.cr <mailto:mspin...@una.ac.cr>
mspinol...@gmail.com
Teléfono: (506) 8706 - 4662
Institutional website: ICOMVIS
<http://www.icomvis.una.ac.cr/index.php/manuel>
Blog sobre Ciencia de Datos: https://mspinola-ciencia-de-datos.netlify.app
--
Marcelino de la Cruz Rot
Coordinador funcional de Biología
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
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 nested grid in R
Dear Manuel, This is R. There is no "it is possible". Only "how" ;-). For example, with spatstat.geom, A <- tess(xgrid=seq(0,80, by=8),ygrid=seq(0,80, by=8)) B <- tess(xgrid=seq(0,80, by=4),ygrid=seq(0,80, by=4)) C <- tess(xgrid=seq(0,80, by=2),ygrid=seq(0,80, by=2)) D <- tess(xgrid=seq(0,80, by=1),ygrid=seq(0,80, by=1)) Cheers, Marcelino El 10/02/2023 a las 0:21, Manuel Spínola escribió: Dear list members, Is it possible to generate a spatial nested grid in R? For example, a grid of several 8km x 8km tiles, and within that grid, I want 4 tiles of 4km x 4km, and in each of those I want 4 tiles of 2km x 2km, and in each of those I want 4 tiles of 1km x 1km. Manuel -- Marcelino de la Cruz Rot Coordinador funcional de Biología Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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 with writing a function that attaches atributes from polygons of a shapefile to the points of a grid
Hi Miguel,
Maybe this way:
xy<-gridcentres(madrid_sqr_m,150,150)
xy.df <- as.data.frame(xy)
xy_sf <- st_as_sf(xy.df, coords=1:2)
xy_crs <- st_sfc(xy_sf, crs=25830)
xy_joined <- st_join(xy_crs, Catfiscal)
join.ppp <- ppp(st_coordinates(xy_joined)[,"X"],
st_coordinates(xy_joined)[,"Y"], window=madrid_sqr_m, marks
=xy_joined$categoria)
Cheers,
Marcelino
El 23/01/2023 a las 20:11, MIGUEL GOMEZ DE ANTONIO escribió:
[No suele recibir correo electr?nico de m...@ccee.ucm.es. Descubra por qu? esto
es importante en https://aka.ms/LearnAboutSenderIdentification ]
Hi all,
I need to create a function that georreference covariates that are not
available at the point level but at a census tract or any other
administrative polygon level (cities, regions, etc.).
My last goal is to create my own function that attaches the values of a
census tract shapefile to the points of the grid in an observation window.
Then I would like to kernel smooth it to obtain a covariate that is
continouos in the observation window that contains the grid points.
I am guessing that my problem is either in creating an object of the class
sf when adding the dataframe to the sfg object or when using the function
st_join () that is supposed to add the attributes of the polygon to the
points located in it.
#I first created the coordinates of the points of the grid in the
observation window "madrid_sqr_m" xy<-gridcentres(madrid_sqr_m,150,150)#
"madrid_sqr_m" is just the boundary of the grid (class owin)
#Then I created a sfg object multipolygon
x=as.numeric(xy$x)
y=as.numeric(xy$y)
XY=data.frame(x,y)
xy=as.matrix(XY)
xy=st_multipoint(xy)#creo el objeto sfg que es un conjunto de puntos
(multypoligon)
xy_crs=st_sfc(xy, crs=25830)#add the proyected reference system
xy_crs
#add the dataframe to have an object of class sf. NOT SURE IF I AM DOING
THIS RIGHT. THE GRID POINTS HAVE NOT ATTRIBUTES SO THE DATAFRAME SHOULD BE
EMPTY, ISN`T?
df <- data.frame(nombre = c("grid"))#
grid_sf=st_sf(df, xy_crs)#
class(grid_sf)# #Add the covariates of a shapefile named "Catfiscal", which
is a census tract polygon, to the points of the grid
join<-st_join(grid_sf,Catfiscal)# PROBLEM MIGHT BE HERE BECAUSE IS NOT
ASIGNING TO EACH GRID POINT THE COVARIATES OF THE POLIGONS
x=join$categoria# "categoria" is the covariate I would like to add to the
points of the grid once I convert the sf object into a "ppp" object of the
Spatstat
join.ppp<-as.ppp(st_coordinates(join), madrid_sqr_m)#"madrid_sqr_m" is just
the boundary of the grid (class owin)
grid_con_cat_f.ppp<-setmarks(join.ppp,x) **
*#Error in `marks<-.ppp`(`*tmp*`, value = value) : # number of points !=
number of marks ** *
density.CATF<-Smooth.ppp(grid_con_cat_f.ppp)
I have been stuck on this several weeks os any advice with be very much
appreciated
Thanks for your help
Miguel
[[alternative HTML version deleted]]
_______________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Marcelino de la Cruz Rot
Coordinador funcional de Biología
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] raster::levels() not working in packaged function.
Maybe including
import(raster)
or
importFrom("raster", "levels")
in the package NAMESPACE would help.
Cheers,
Marcelino
El 31/10/2020 a las 13:24, nevil amos escribió:
Apologies, I cannot see how to make a rero for this issue.
I have a function that uses levels(r) tor return the RAT of a raster "r"
when the function is sourced from a script
source(".\R\function.r")
it works fine.
when the function is built into a package and sourced from there
library(mypackage) using the same script file to make the package
levels(r)[[1]]
the same line throws an error, as levels(myraster returns NULL
If I modify the script to include the raster namespace:
raster::levels(r)[[1]]
Then I get the error
Error: 'levels<-' is not an exported object from 'namespace:raster'
I have also tried just using levels(r) and putting raster as a depends
rather than an import in the DESCRIPTION file for the package, this does
not solve the error.
Any suggestions on how to overcome the problem?
many thanks
Nevil Amos
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] unsubscribing R mailing list
It is very easy: Go to the end of this web page: https://stat.ethz.ch/mailman/listinfo/r-sig-geo and unsubscribe yourself from the list. Cheers, Marcelino El 21/09/2020 a las 18:58, Adeela Munawar escribió: How to unsubscribe this mailing list? [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo . -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Filtering a set of points in a "ppp" object by distance using marks
Sorry, I meant to say
subset(insects.ppp, marks=="termiNests" & ddd[,"antNests"] >20)
El 16/09/2020 a las 9:18, Marcelino de la Cruz Rot escribió:
Hi Alexandre,
may be this?
ddd <- nndist(insects.ppp, by=factor(insects.ppp$marks))
subset(insects.ppp, marks=="termiNests" & ddd[,"termiNests"] >20)
Cheers,
Marcelino
El 15/09/2020 a las 22:52, ASANTOS via R-sig-Geo escribió:
Dear R-Sig-Geo Members,
I'd like to find any way to filtering a set of points in a "ppp"
object by minimum distance just only between different marks. In my
example:
#Package
library(spatstat)
#Point process example - ants
data(ants)
ants.ppp<-ppp(x=ants$x,y=ants$y,marks=rep("antNests",length(ants$x)),window=Window(ants))
# Create a artificial point pattern - termites
termites <- rpoispp(0.0005, win=Window(ants))
termites.ppp<-ppp(x=termites$x,y=termites$y,marks=rep("termiNests",length(termites$x)),window=Window(ants))
#Join ants.ppp and termites.ppp
insects.ppp<-superimpose(ants.ppp,termites.ppp)
#If I try to use subset function:
subset(insects.ppp, pairdist(insects.ppp) > 20 & marks=="termiNests")
#Marked planar point pattern: 223 points #marks are of storage type
�character� #window: polygonal boundary #enclosing rectangle: [-25, 803]
x [-49, 717] units (one unit = 0.5 feet) #Warning message: #In ppp(X[,
1], X[, 2], window = win, marks = marx, check = check) : # 70751 out of
70974 points had NA or NaN coordinate values, and were discarded
Not the desirable result yet, because I'd like to calculate just only
the > 20 "termiNests" to "antNests" marks and not the "termiNests"
with "termiNests" too.
Please any ideas?
Thanks in advanced,
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Filtering a set of points in a "ppp" object by distance using marks
Hi Alexandre,
may be this?
ddd <- nndist(insects.ppp, by=factor(insects.ppp$marks))
subset(insects.ppp, marks=="termiNests" & ddd[,"termiNests"] >20)
Cheers,
Marcelino
El 15/09/2020 a las 22:52, ASANTOS via R-sig-Geo escribió:
Dear R-Sig-Geo Members,
I'd like to find any way to filtering a set of points in a "ppp" object by
minimum distance just only between different marks. In my example:
#Package
library(spatstat)
#Point process example - ants
data(ants)
ants.ppp<-ppp(x=ants$x,y=ants$y,marks=rep("antNests",length(ants$x)),window=Window(ants))
# Create a artificial point pattern - termites
termites <- rpoispp(0.0005, win=Window(ants))
termites.ppp<-ppp(x=termites$x,y=termites$y,marks=rep("termiNests",length(termites$x)),window=Window(ants))
#Join ants.ppp and termites.ppp
insects.ppp<-superimpose(ants.ppp,termites.ppp)
#If I try to use subset function:
subset(insects.ppp, pairdist(insects.ppp) > 20 & marks=="termiNests")
#Marked planar point pattern: 223 points #marks are of storage type
�character� #window: polygonal boundary #enclosing rectangle: [-25, 803]
x [-49, 717] units (one unit = 0.5 feet) #Warning message: #In ppp(X[,
1], X[, 2], window = win, marks = marx, check = check) : # 70751 out of
70974 points had NA or NaN coordinate values, and were discarded
Not the desirable result yet, because I'd like to calculate just only the > 20 "termiNests" to
"antNests" marks and not the "termiNests" with "termiNests" too.
Please any ideas?
Thanks in advanced,
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Simple Ripley's CRS test for market point patters
Dear Alexandre, I think that the solution to your first problem would be something like this: #Area 0-25 env.formi1<-envelope(syn.ppp[syn.ppp$marks>=0 & syn.ppp$marks<25], nsim=99,fun=Kest) plot(env.formi1,lwd=list(3,1,1,1), main="") Then, for your second problem, you should clarify what do you understand by "market intensity observed pattern". In fact, rmpoispp() simulates random multitype ppp's, i.e., it does not estimate parameters of quantitatively marked point processes. Maybe be you would consider interesting Smooth.ppp(syn.ppp) and envelope(syn.ppp, nsim=99,fun=markcorr, simulate=expression(rlabel(syn.ppp))) Cheers, Marcelino El 22/07/2019 a las 23:36, Alexandre Santos via R-sig-Geo escribió: > Dear R-Sig-Geo Members, I"ve like to find any simple way for apply CRS > test for market point patters, for this I try to create a script below: > #Packages require(spatstat)require(sp) > > # Create some points that represents ant nests in > UTMxp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,371159.291,371158.552,370978.252,371120.03,371116.993) > yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926) > # Then I create the size of each nest - my covariate used as marked > processarea<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,105,133,135,23.92,190,12.9,15.2,192.78,104,255,24) > > # Make a countour - only as exerciseW <- convexhull.xy(xp,yp) > #Create a ppp > objectp_xy<-cbind(xp,yp)syn.ppp<-ppp(x=coordinates(p_xy)[,1],y=coordinates(p_xy)[,2],window=W, > marks=area)syn.ppp <- as.ppp(syn.ppp)plot(syn.ppp, main=" ") > First, I've like to study CSR of market point process (my hypothesis is that > different size have a same spatial distribution) when area >= 0, area < 25 > and area >=25, area < 55, for this I make: > # Area 0-25env.formi1<-envelope(syn.ppp,nsim=99,fun=Kest, area >= 0, area < > 25)plot(env.formi1,lwd=list(3,1,1,1), main="") > # Area 25-55env.formi2<-envelope(syn.ppp,nsim=99,fun=Kest, area >=25, area < > 55)plot(env.formi2,lwd=list(3,1,1,1), main="") > My first problem is that I have the same plot in both conditions and I don't > know why. > Second, if I try to estimate the market intensity observed pattern > est.int <- density(syn.ppp)est_xy <- rmpoispp(est.int)plot(est_xy, main=" ") > My output is only my points position without marked area in my ppp object > created. > My question is what is the problem with this Ripley's reduced second moment > function approach? There are any way for study my point process when the area > is a covariate of my point process? > Thanks in advanced > Alexandre > > > -- > ==Alexandre > dos SantosProteção Florestal IFMT - Instituto Federal de Educação, Ciência e > Tecnologia de Mato GrossoCampus CáceresCaixa Postal 244Avenida dos Ramires, > s/nBairro: Distrito Industrial Cáceres - MT CEP: > 78.200-000Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 > (FIXO)e-mails:alexandresanto...@yahoo.com.br > alexandre.san...@cas.ifmt.edu.br Lattes: > http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 > - ResearcherID: A-5790-2016Researchgate: > www.researchgate.net/profile/Alexandre_Santos10 > LinkedIn: > br.linkedin.com/in/alexandre-dos-santos-87961635Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ > ================== > > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > . > -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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 with enveloped test in spatstat
Hi David, This is very clearly explained in page 393 of Baddeley et al. 2015. Basically, the MAD test is affected by transformations of the summary function, with the L-function providing more powerful tests because of its stabilization of the variance. If you are interested in point pattern analysis I would recommend you to get a copy of this nice book. Cheers, Marcelino Adrian Baddeley, Ege Rubak, Rolf Turner (2015). Spatial Point Patterns: Methodology and Applications with R. London: Chapman and Hall/CRC Press. http://www.crcpress.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/9781482210200/ El 21/04/2018 a las 18:05, David Unwin escribió: Can any *spatstat* user explain to me why the two p-values obtained below for an envelope test against CSR are so different? data(swedishpines) d<-swedishpines plot(d) mad.test(d,Kest,nsim=999,verbose=F) Maximum absolute deviation test of CSR Monte Carlo test based on 999 simulations Summary function: K(r) Reference function: theoretical Alternative: two.sided Interval of distance values: [0, 24] units (one unit = 0.1 metres) Test statistic: Maximum absolute deviation Deviation = observed minus theoretical data: d mad = 150.69, rank = 216, p-value = *0.216* mad.test(d,*Lest*,nsim=999,verbose=F) Maximum absolute deviation test of CSR Monte Carlo test based on 999 simulations Summary function: L(r) Reference function: theoretical Alternative: two.sided Interval of distance values: [0, 24] units (one unit = 0.1 metres) Test statistic: Maximum absolute deviation Deviation = observed minus theoretical data: d mad = 2.9921, rank = 9, p-value = *0.009* *!!!* These data are dispersed relative to CSR: kl<-envelope(d,Kest,nsim=999,correction="border") plot(kl) Dave Unwin [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo . -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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] How to label a polygon window (spatstat package)
Maybe you meant something like this:
# example of several polygons (in a list)
x <- runif(20)
y <- runif(20)
polygons<- list()
for (i in 1:20) polygons[[i]]<- disc(radius=0.05, centre=c(x[i], y[i]))
plot(owin())
lapply(polygons, plot, add=T)
# label each polygon
labs<-letters
for (i in 1: length(polygons)) {xy<-centroid.owin(polygons[[i]]);
text(xy$x,xy$y, labs[i])}
Cheers,
Marcelino
El 10/02/2018 a las 6:43, Mohammad Tanvir Ahamed via R-sig-Geo escribió:
Thanks.
when there is multiple polygon , it a problem .
looking for something more .
Regards.
Tanvir Ahamed
Stockholm, Sweden | mashra...@yahoo.com
On Saturday, February 10, 2018, 6:35:59 AM GMT+1, Michael Sumner
<mdsum...@gmail.com> wrote:
Try
text(0.5, 0.5, label = "?text")
On Sat, 10 Feb 2018, 16:22 Mohammad Tanvir Ahamed via R-help,
<r-h...@r-project.org> wrote:
Hi,
I want to label a polygon (circle or polygon) inside.
As for example code
library(spatstat)
x <- runif(20)
y <- runif(20)
X <- ppp(x, y, window=disc(0.7))
plot(X)
Now I want to label that circle inside . Can some one please help me ?
Thanks.
Regards.
Tanvir Ahamed
Stockholm, Sweden | mashra...@yahoo.com
__
r-h...@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Splitting spatial lines at self intersections (line crossings)
You should load before library(maptools) and the coercion is done. Cheers Marcelino El 26/01/2018 a las 13:12, Michael Sumner escribió: Ege: I don't believe that this code is enough to get coercion between sp and spatstat: pts <- cbind(c(120:123,121,125),c(100,100,104,102,99,98)) pt2 <- cbind(c(124,124,123,118,124,125),c(100,97,100,104,106,110)) library(sp) library(spatstat) L1 <- SpatialLines(list(Lines(list(Line(pts)), "X"))) L2 <- SpatialLines(list(Lines(list(Line(pt2)), "X"))) l1 <- as.psp(L1) l2 <- as.psp(L2) int <- crossing.psp(l1, l2) Are we relying on some other code, or perhaps a specific version of spatstat or its family? (I know ways to convert from Spatial to psp, but I was taken by this apparent lack of infrastructure - and I'm very keen on seeing more bridges between these rich worlds). Cheers, Mike On Fri, 26 Jan 2018 at 17:48 Ege Rubak <ru...@math.aau.dk> wrote: The `psp` class in spatstat consists of individual line segments and `selfcrossing.psp` checks whether each individual line intersects one of the other lines, which happens at all the points in your plot. If instead you treat each of the two line sequences as a `psp` you can check find the crossings of the two `psp` objects. I.e., continuing your code (with `spatstat` and `maptools` loaded): L1 <- SpatialLines(list(Lines(list(Line(pts)), "X"))) L2 <- SpatialLines(list(Lines(list(Line(pt2)), "X"))) l1 <- as.psp(L1) l2 <- as.psp(L2) int <- crossing.psp(l1, l2) This gives you the four intersections between the lines. I don't know of a simple way to do the next task in `spatstat`. A useful function if you are going to try to write some code is `test.crossing.psp` which gives you a logical matrix indicating which segments cross. Cheers, Ege On 01/25/2018 11:26 PM, Glenn Stauffer wrote: I have a Spatial Lines object I would like to split at every point where the line self-intersects (crosses or touches itself), or anywhere lines touch each other, if there are multiple lines. I've tried using spatstat to convert the SpatialLines to a psp object and then use the selfcrossing.psp function to find the intersection points (see example below). But that seems to identify all the nodes, not just the points where the line crosses. Also, even if it identified only the crossings (which is what I want), I am not sure how to perform the next step and split the lines at those points. So, 1) I need to identify the crossings/touches, and 2) split the lines at those points. Essentially, what I am looking for is an R equivalent to the planarize function in ArgGIS. Is there a relatively easy way to do this in R? Thanks, Glenn Here is an example of what I have tried. pts <- cbind(c(120:123,121,125),c(100,100,104,102,99,98)) pt2 <- cbind(c(124,124,123,118,124,125),c(100,97,100,104,106,110)) projstr <- "+init=epsg:3071" # make everything in meters L <- SpatialLines(list(Lines(list(Line(pts),Line(pt2)),"X")),proj4string = CRS(projstr)) plot(L) PSP <- as.psp.SpatialLines(L) int <- selfcrossing.psp(PSP) plot(int,add=TRUE,col="red") # identifies more than just the crossings [[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 -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Autologistic regression in R
Hi Mingke Li, Dorman et al. (2007) review different methods to account for spatial dependence in linear models, including the use of autocovariables, and provide R code for them. HTH, Marcelino Dormann, et al. 2007. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. – Ecography 30: 609–628. http://www.ecography.org/appendix/e5171 El 15/11/2017 a las 2:46, Mingke Li escribió: Hi, I am new to autologistic regression and R. I do have questions when starting a project in which I believe autologistic regression (spdep package) is needed. I have a point layer whose attribute table stores the values of the dependent variable (population of a kind of insect), all the independent variables (environmental factors), and the associated latitude and longitude. I hope to to fit an autologistic model to analyze which factors or combinations of factors have effects on the presence/absence of the insect (1 or 0). I found other papers which applied autologistic regression in their study almost used a grid system and defined their window sizes. So, my question is do I have to convert my point layer into a grid system if I want to do this analysis with R? Also, what should I consider when I generate the grid system? How to determine a proper cell size? How about the searching window (neighbourhood) size? Many Thanks. Erin [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo . -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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 clusters with equal number of objects (2)
Hello Maurizio, If they are cells in a grid, they are already as much clustered as they could be ;-p. Please, be more specific. Cheers, Marcelino El 31/10/2017 a las 11:57, Maurizio Marchi escribió: Hello everybody, is there an easy way to cluster spatial objects (polygons in my case, a grid of 500mx500m cells) with an equal number of objects per cluster? I tried to attach an example of what I mean but I think the message was never delivered. Regards -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] WorldClim Version2
I'm afraid it don't.
At least in version 2.5.8, as it appears from this line in .wordclim():
theurl <-
paste("http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/;,
zip, sep = "")
Marcelino
El 28/09/2017 a las 4:31, Bacou, Melanie escribió:
I would assume `raster::getData("worldclim", ...)` uses WorldClim
version 2?
--Mel.
On 09/25/2017 05:14 PM, Andy Bunn wrote:
Great. Thanks Marcelino. I have all the data dowloaded. I'm surprised that
there isn't a dedicated package for worldclim yet ‹ making one that ties
in with packages like BIEN would be great. Perhaps I'll work on it someday.
-A
On 9/25/17, 11:32 AM, "Marcelino de la Cruz Rot"
<marcelino.delac...@urjc.es> wrote:
Hello Andy,
I've been using raster to play with WordlClim 2.0 and it works fine. For
example:
# Download 2.5' arc P data of worldclim 2.0
# The file is 68.5 MB so it takes some time to download
P_url<-
"http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_2.5m_prec.zi
p"
download.file(P_url, destfile="wc2.0_2.5m_prec.zip")
unzip("wc2.0_2.5m_prec.zip")
# Generate a stack object with all monthly precipitation layers
P_filenames<- paste("wc2.0_2.5m_prec_", c(paste(0,1:9, sep=""), 11,12),
".tif",sep="")
prec2.5<- stack(P_filenames)
# Compute annual precipitation as a raster layer
P2.5<- sum(prec2.5)
# etc.
Cheers,
Marcelino
El 25/09/2017 a las 19:42, Andy Bunn escribió:
Hello all, is anybody aware of a R package that gracefully works with
version 2 of the WorldClim data?
http://worldclim.org/version2
I have all the data extracted locally but I'm wondering if there is
something akin to getData in raster that will work with version 2 and
not
1.4 as (the brilliant) raster package currently does.
Many thanks, Andy
___________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] WorldClim Version2
Hello Andy,
I've been using raster to play with WordlClim 2.0 and it works fine. For
example:
# Download 2.5' arc P data of worldclim 2.0
# The file is 68.5 MB so it takes some time to download
P_url<-
"http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_2.5m_prec.zip;
download.file(P_url, destfile="wc2.0_2.5m_prec.zip")
unzip("wc2.0_2.5m_prec.zip")
# Generate a stack object with all monthly precipitation layers
P_filenames<- paste("wc2.0_2.5m_prec_", c(paste(0,1:9, sep=""), 11,12),
".tif",sep="")
prec2.5<- stack(P_filenames)
# Compute annual precipitation as a raster layer
P2.5<- sum(prec2.5)
# etc.
Cheers,
Marcelino
El 25/09/2017 a las 19:42, Andy Bunn escribió:
Hello all, is anybody aware of a R package that gracefully works with
version 2 of the WorldClim data?
http://worldclim.org/version2
I have all the data extracted locally but I'm wondering if there is
something akin to getData in raster that will work with version 2 and not
1.4 as (the brilliant) raster package currently does.
Many thanks, Andy
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] quadracount on multitype points
Hi Guy,
I only like real pies :-) but, what about something like this:
library(spatstat)
data(lansing)
qs<- quadratcount(split(lansing))
dqs <- dim(qs[[1]])
nr <- dqs[1]
nc <- dqs[2]
le <- length(qs[[1]])
m<- matrix( 1:le, nr=nr, nc=nc, byrow=T)
layout(m, heights=rep(1,nr), widths=rep(1,nc))
#layout.show(25)
par(mar=c(0,0,0,0))
for (i in 1:nc){
for (j in 1:nr){
pie(sapply(qs, function(x) x[i,j]), labels=NA )
box()
}
}
Cheers,
Marcelino
El 13/09/2017 a las 10:21, Ege Rubak escribió:
Hi Guy,
Maybe your explorative analysis could also benefit from `relrisk` in
case you don't know that function?
It basically gives you a list of smooth images (one for each type) of
the probability that a hypothetical sample at that location is of a
given type. E.g. for the built-in multitype point pattern dataset
`lansing` you would do:
library(spatstat)
prob <- relrisk(lansing, diggle=TRUE)
plot(prob)
This example and more (e.g. a division of the sampling area into
subsets of dominant types) can be found in Chapter 14 of our "spatstat
book" if you need more details (sorry about the shameless self
promotion, but I don't know a better source for this :-)).
Kind regards,
Ege
On 09/13/2017 01:08 AM, Guy Bayegnak wrote:
Thanks a lot for your response and suggestion Rolf. Yes, by
"quadrant" I mean the little sub-windows. My problem is the following:
We have collected thousands of groundwater samples across a vast
area, and analysed them. Based on the analysis we are able to assign
a "type" to each water sample. When plotted, there seems to be a
spatial trend in water type. But a given area may have more than one
water type, usually with a dominant type (most frequently
occurring). What I am trying to do is identify the dominant type for
each sub-region /sub-windows but show the count side by side, for
example:
x
y [0,0.801) [0.801,1.6]
[0.5,1] Off = 36 Off = 6
On = 3 On = 39
[0,0.5) Off = 4 Off = 36
On = 42 On = 6
I think I can achieve what I am looking for with your suggestion.
Once I get the table list, I will copy the numbers side by side
manually.
Sincerely,
Guy
-Original Message-
From: Rolf Turner [mailto:r.tur...@auckland.ac.nz]
Sent: September 12, 2017 3:45 PM
To: Guy Bayegnak <guy.bayeg...@gov.ab.ca>
Cc: r-sig-geo@r-project.org; adrian.badde...@curtin.edu.au; Ege Rubak
<ru...@math.aau.dk>
Subject: Re: [R-sig-Geo] quadracount on multitype points
On 13/09/17 02:11, Guy Bayegnak wrote:
Dear all,
I am working on a multitype point pattern, and I'd like to estimate
how many of each type of point occurs into each quadrant. I know it is
possible to use the quandracount on split marks as follows using
spatstats: quadratcount(split(marks)). But the result produces as many
windows as they are marks. I am wondering is there is a way to know
many occurrence of each type there is per quadrant and to plot it in a
single grid.
Thanks.
You really should start by mentioning that you are dealing with the
spatstat package.
It's not clear to me what you are after. A minimal reproducible
example would be helpful. I presume that by "quadrant" you mean one
of the four equal sub-windows formed by bisecting your (rectangular)
window vertically and horizontally.
If my presumption is correct then perhaps
lapply(split(X),quadratcount,nx=2)
(where "X" is your point pattern) does what you want. E.g.:
lapply(split(amacrine),quadratcount,nx=2)
$off
x
y [0,0.801) [0.801,1.6]
[0.5,1] 36 36
[0,0.5) 34 36
$on
x
y [0,0.801) [0.801,1.6]
[0.5,1] 35 39
[0,0.5) 42 36
Is this something like what you wish to achieve?
cheers,
Rolf Turner
--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
This email and any files transmitted with it are confi...{{dropped:3}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] quadracount on multitype points
Hi Guy,
You can know how many points of each type are in each cell with sapply,
i.e.:
>data (lansing)
>qs<- quadratcount(split(lansing))
#trees in the fist row, first column:
> sapply(qs, function(x) x[1,1])
blackoak hickory maple misc redoak whiteoak
13 66 0 0 4 23
#trees in the fist row, second column:
> sapply(qs, function(x) x[1,2])
blackoak hickory maple misc redoak whiteoak
16 38 3 0 9 16
# ETC. On the other hand I don't know if it is a good idea to have all
this numbers printed in just one cell or how to print them without the
result being a mess...
Cheers,
Marcelino.
El 12/09/2017 a las 16:11, Guy Bayegnak escribió:
Dear all,
I am working on a multitype point pattern, and I'd like to estimate how many of
each type of point occurs into each quadrant. I know it is possible to use the
quandracount on split marks as follows using spatstats:
quadratcount(split(marks)). But the result produces as many windows as they are
marks. I am wondering is there is a way to know many occurrence of each type
there is per quadrant and to plot it in a single grid.
Thanks,
Guy
This email and any files transmitted with it are confide...{{dropped:10}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
.
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Wavelet analysis of anisotropy in point patterns
Hi Domenico, AFAIK, you should build your own functions to get the same analysis provided by the PASSaGE software. If you have some elementary R programming skills, following the description of the analysis proposed by Rosenberg(2004) shouldn't be much complicated. Otherwise...I don't know. Good luck! Marcelino El 19/07/2017 a las 12:12, Domenico Giusti escribió: Thank you Marcelino, I already had a look at those functions and, although they do the job, I was wondering how to analyze anisotropy in point pattern using wavelet analysis. I tried to use the WaveletComp package, but with no success. It looks to me that packages for wavelet analysis work fine with temporal data, but I couldn't find the way to use it with spatial point patterns. How could I get in R the same analysis provided by the PASSaGE software? Thank you! On 07/18/2017 06:17 PM, Marcelino de la Cruz Rot wrote: Hi Domenico, In spatstat you have several functions to detect anysotropy in point patterns: - the sector K-function, Ksector(); - the pair orientation distribution, pairorient(); - the anisotropic pair correlation function, Kmeasure(). See the help pages of these functios or, even better, consult pages 236:242 in Adrian Baddeley, Ege Rubak, Rolf Turner (2015). Spatial Point Patterns: Methodology and Applications with R. London: Chapman and Hall/CRC Press, 2015. http://www.crcpress.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/9781482210200/ Cheers, Marcelino El 18/07/2017 a las 16:30, Domenico Giusti escribió: Hi all, I'm running wavelet analysis of anisotropy in point patterns using the PASSaGE software by Rosenberg et al. Rosenberg, M.S., and C.D. Anderson (2011) PASSaGE: Pattern Analysis, Spatial Statistics and Geographic Exegesis. Version 2. Methods in Ecology & Evolution 2(3):229-232. http://www.passagesoftware.net Rosenberg, M.S. 2004. Wavelet analysis for detecting anisotropy in point patterns. Journal of Vegetation Science 15:277-284. DOI: 10./j.1654-1103.2004.tb02262.x I wonder how I can do the same analysis using R. Thanks, Dome ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Wavelet analysis of anisotropy in point patterns
Hi Domenico, In spatstat you have several functions to detect anysotropy in point patterns: - the sector K-function, Ksector(); - the pair orientation distribution, pairorient(); - the anisotropic pair correlation function, Kmeasure(). See the help pages of these functios or, even better, consult pages 236:242 in Adrian Baddeley, Ege Rubak, Rolf Turner (2015). Spatial Point Patterns: Methodology and Applications with R. London: Chapman and Hall/CRC Press, 2015. http://www.crcpress.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/9781482210200/ Cheers, Marcelino El 18/07/2017 a las 16:30, Domenico Giusti escribió: Hi all, I'm running wavelet analysis of anisotropy in point patterns using the PASSaGE software by Rosenberg et al. Rosenberg, M.S., and C.D. Anderson (2011) PASSaGE: Pattern Analysis, Spatial Statistics and Geographic Exegesis. Version 2. Methods in Ecology & Evolution 2(3):229-232. http://www.passagesoftware.net Rosenberg, M.S. 2004. Wavelet analysis for detecting anisotropy in point patterns. Journal of Vegetation Science 15:277-284. DOI: 10./j.1654-1103.2004.tb02262.x I wonder how I can do the same analysis using R. Thanks, Dome -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Test for difference between two pint patterns
Hi Paolo, see the help page of Kcross(), pcfcross() or any other crossed function in spatstat (and also that of envelope() for the tests). Best, Marcelino El 17/07/2017 a las 10:32, Paolo Piras escribió: Hi folks, I'm new in point pattern analysis; thus my question could appear trivial for experts; Given two ppp objects (spatstat package) [with different number of points] I'm looking for a way to assess the statistical difference in the distributions of points. studpermu.test() seems doing the job but for groups of point patterns while I have only two distributions. Kest() could be used to look for deviation from Poisson distribution but I should compare the two patterns. I think a hand-made permutation test where points are randomly re-assigned between the two distribution and K stat is computed at any run (and compared with the observed ones) could be ok but I'm not sure however. Thanks in advance for any advice Paolo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo . -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] [FORGED] Re: Making an hexagonal grid using spsample
My apologies for such a pair (or more) of embarrassing mistakes! I should read a bit more these days... Marcelino El 04/03/2017 a las 2:31, Rolf Turner escribió: On 04/03/17 08:38, Marcelino de la Cruz Rot wrote: Hi Manuel, I do answer to the question "How can I make a spatial grid of 1 ha (or other size) in R?" You can use function hextess in spatstat library(spatstat) # some arbitrary area, with coordinates in hectometres W <- Window(chorley) # As Rolf said, hexagons of 1ha should have side of 402.0673 metres, so, in hectometres: s <- 4.020673 plot(hextess(W, s)) plot(hexgrid(W, s), add=TRUE) Marcelino, Actually I said hexagons of area *42* ha should have side length equal to 402.0673 metres. Moreover the Chorley data set has units of *kilometres* not hectometres, so that should be s <- 0.4020673. Or, to avoid just a touch of round-off error, s <- sqrt(2*0.42)/3^0.75. Note that if you then do xxx <- hextess(W,s,trim=FALSE) unique(sapply(tiles(xxx),area.owin)) you get 0.42 --- i.e. 0.42 square kilometres, or 42 hectares. cheers, Rolf -- Marcelino de la Cruz Rot Depto. de Biología y Geología Física y Química Inorgánica Universidad Rey Juan Carlos Móstoles España ___ 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 an hexagonal grid using spsample
Hi Manuel,
I do answer to the question "How can I make a spatial grid of 1 ha (or
other size) in R?"
You can use function hextess in spatstat
library(spatstat)
# some arbitrary area, with coordinates in hectometres
W <- Window(chorley)
# As Rolf said, hexagons of 1ha should have side of 402.0673 metres, so,
in hectometres:
s <- 4.020673
plot(hextess(W, s))
plot(hexgrid(W, s), add=TRUE)
Cheers,
Marcelino
El 03/03/2017 a las 20:03, Manuel Spínola escribió:
Dear list members,
I am trying to make an spatial hexagonal grid in R using the spsample
function. I would like hexagons of 1 ha = 1 sq meters.
First I used a cell size of 62.04 m as the side of the hexagon.
data(meuse.grid)
gridded(meuse.grid) = ~x + y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04)
xx <- HexPoints2SpatialPolygons(xx)
rgeos::gArea(xx, byid=TRUE)
It gave me an area per hexagon = .299
Second I used a cell size of 62.04 * 2 = 124.08 m as the long diagonal of
the hexagon
data(meuse.grid)
gridded(meuse.grid) = ~x + y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04*2)
xx <- HexPoints2SpatialPolygons(xx)
rgeos::gArea(xx, byid=TRUE)
It gave me an area per hexagon = 1.19
Both results are wrong, because I was expecting close to 1.
How can I make a spatial grid of 1 ha (or other size) in R?
Manuel
--
Marcelino de la Cruz Rot
Depto. de Biología y Geología
Física y Química Inorgánica
Universidad Rey Juan Carlos
Móstoles España
___
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 point density raster OR other ways to thin point data
Hi, you can combine density.ppp and rthin to thin inversely to the density of points. See this example # example of a heterogeneous point pattern library(spatstat) data(lansing) maple-lansing[lansing$marks==maple] plot(maple) maple.dens- density(maple) plot(maple.dens) # get the probability of retention as inversely proportional to the density maple.p - eval.im(1-(maple.dens/max(maple.dens))) plot(maple.p) # Thin your pp plot(rthin(maple, maple.p)) Sigma us calculted by default by density.ppp but you can set your aown sigma also (see the help page of density.ppp) HTH, Marcelino El 2013-02-24 19:38, Julie Lee-Yaw escribió: Hi I have a set of sampling locations with a strong bias towards one area within my study region (e.g. many more points in the south than in the north). I want to thin sites in this high density region by randomly discarding sites until the density of sampling is similar to other areas of my study region. The approach I've come up with so far is something like: 1) calculate a kernel density or point density raster from the points 2) take the average density value of all points above (x) and below (y) a user-defined cutoff 3) for all points with an initial density value a defined cutoff, randomly discard point one at a time, create a new density raster, recalculate the average density (x) of these points...repeat until x = y To make this work, I'm trying to figure out how to produce a kernel density raster or point density raster from point data in R (akin to the kernel density and point density tools in ArcGIS). Can anyone point me in the right direction? I have seen the density.ppp function but I am unsure of how to calculate sigma or whether this function even what I need to meet my ultimate goal. Any alternative solutions for thinning point data are also appreciated. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- MARCELINO DE LA CRUZ ROT Universidad Politecnica de Madrid ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Contour plot for point process
I don't know if I uuderstand you but maybe this is what you are asking for: xypois.d - density(xypois) plot(xypois.d, col=grey((0:100)/100)) contour(xypois.d, add=T) points(xypois, pch=19, cex=0.7) See the help page of density.ppp to customize the interpolation HTH, MArcelino El 2013-02-05 12:32, ASANTOS escribió: Dear list members, I'm trying to create a contour plot in gray scale and convert in pixel format for represent a interpolation of point process, my short script: require(spatstat) A-316## Size of square area # Create points xypois=rMatClust(0.005,7,8,owin(c(0,A),c(0,A))) xd=xypois$x yd=xypois$y plot(xd,yd,xlim=c(0,A),ylim=c(0,A),pch=19,cex=0.5)#Draw # I ask, if there are any function for make this? Thanks -- Alexandre dos Santos Forest entomology IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- MARCELINO DE LA CRUZ ROT Universidad Politecnica de Madrid ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
