[R-sig-Geo] Clustering spatial points in a SpatialDataFrame object
Dear r-sig-geo Members, I'd like to calculate the centroid of each cluster (considering points with < 10m distance to the same cluster) of spatial points using mean operation (spdplyr package) for coordinates and another operation (sum) for the attribute(area) without success. In my example: #Packages library(sp) library(maptools) library(spdplyr) library(cluster) # Small sample (40 points) small.sample<-read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/sample_points.csv";) #Convert to spatial object xy.small.sample <- small.sample[,c(1,2)] spdf.small.sample <- SpatialPointsDataFrame(coords = xy.small.sample, data = small.sample, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) # Transform to UTM utm.wgs.84 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs" small.sample.utm <- spTransform(spdf.small.sample, utm.wgs.84) # #Convert each cluster to one point Rc=10 #Maximum distance between points is 10 meters small.sample.utm.dm<-spDists(small.sample.utm) # Distance matrix clusters <- as.hclust(agnes(small.sample.utm.dm, diss = T)) #Hierarchical Clustering small.sample.utm@data$class<- cutree(clusters, h=Rc) # Cut into groups each 10 meters # Average of x and y coordinates and area using spdplyr package small.sample.utm.classification<- small.sample.utm %>% group_by (class) %>% summarise (area_clu=mean(area),x_clu=mean(coords[,1]),y_clu=mean(coords[,2])) #Error: Problem with `summarise()` input `x_clu`. #i Input `x_clu` is `mean(coords[, 1])`. #i The error occurred in group 1: class = 1. Here I believe that this error is a cause of the spatial data frame attributes have individual names and coordinates not. If I try to use something like coords[,1], it doesn't work!! My goal is: # Original points representation plot(small.sample.utm, pch=16) # Center of the centroids representation points(small.sample.utm.classification$x_clu,small.sample.utm.classification$y_clu, col="red") # Labelling the area text(small.sample.utm.classification$x_clu ~ small.sample.utm.classification$x_clu, labels=small.sample.utm.classification$area_clu, cex=0.9, font=2) Thanks in advanced, Alexandre ___ 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 pixels neighborhood in a raster [solved]
.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 09/11/2020 19:41, Ben Tupper escreveu: Hi, Your example isn't fully reproducible. Or maybe it just has some little errors. I can't get past this spot... ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) coordinates(ants.res) <- ~ x + y You have requested that raster::rasterToPoints() return a spatial points object, so there shouldn't be a need for sp::coordinates() to cast to a spatial object. Also, I think you are invoking "antscounts" and "antscount" as the same variable. Maybe? It's not clear to me. Finally, when you call raster::adjacent() you haven't defined the value for the cells argument when computing e1 and e2. What do you want them to be? I wish I could be more helpful. Cheers, Ben On Mon, Nov 9, 2020 at 12:46 PM ASANTOS via R-sig-Geo < r-sig-geo@r-project.org> wrote: Dear r-sig-geo Members, I'd like to find any way to create 1 (total 9 pixels) and 2 pixels (total 25 pixels) surrounding the neighborhood of each pixel (ant) in my plot image (antscount). In my example: #Packages library(spatstat) library(raster) #Selection of ants data set data(ants) geo.form<-cbind(x=ants$x,y=ants$y) #Definition of raster resolution - 10 units ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) coordinates(ants.res) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(ants.res) <- TRUE #Rasterize antscount<- rasterize(geo.form, raster(ants.res), fun='count', background=0) values(antscount)[values(antscount) > 0] = 1 #Vizualize plot(antscounts) Now, the selection of neighborhood pixels sounds easy, something like: # For 1 pixel neighborhood neigh1 <- matrix(1L, nrow=3, ncol=3) e1<-adjacent(antscounts, cells , directions=neigh1, pairs=FALSE) ng_coords1 <- xyFromCell(antscounts, e1) # For 2 pixel neighborhood neigh2 <- matrix(1L, nrow=5, ncol=5) e2<-adjacent(antscounts, cells , directions=neigh2, pairs=FALSE) ng_coords5 <- xyFromCell(antscounts, e2) But for the combination of all the information (0 and 1 new pixel values) and the new raster representation (antscounts) I don't have success. Please, any ideas? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ 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] Create pixels neighborhood in a raster
Dear Ben, Sorry about of lot of mistakes, you are sure some part of code are confusing and/or with errors. I'll try again with some corrections: First, I have the rasterization process of ants presence (0 is absence) in a 10x10 units raster: #Packages library(spatstat) library(raster) #Selection of ants data set data(ants) geo.form<-cbind(x=ants$x,y=ants$y) #Definition of raster resolution - 10 units ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) # coerce to SpatialPixelsDataFrame gridded(ants.res) <- TRUE #Rasterize antscount<- rasterize(geo.form, raster(ants.res), fun='count', background=0) values(antscount)[values(antscount) > 0] = 1 #Vizualize plot(antscount) Now, I'd like to create 1 pixels around each pixel of ant presence (Total of 9 pixels in each ant presence in antscount raster): # For 1 pixel neighborhood neigh1 <- matrix(1L, nrow=3, ncol=3); neigh1[2,2] <- 0L ants1<-which(values(antscount)> 0) cells<- xyFromCell(antscount, ants1) e1<-adjacent(antscount, cells, directions=neigh1, pairs=FALSE) ng_coords1 <- xyFromCell(antscount, e1) points(ng_coords1, col="red") I need that's this new create pixels has 1 as value too. #Rasterize for 1 pixel neighborhood antscount.9<- rasterize(ng_coords1 , raster(ants.res), fun='count', background=0) plot(antscount.9) The problem is my ng_coords1 coordinates is wrong and just only in the top of the antscount raster, despite xyFromCell(antscount, ants1) condition. My goal is a new ant presence raster with 8 pixels surrounding the neigourhood of each pixel (ant) in the original antscount raster. Any ideas? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 09/11/2020 19:41, Ben Tupper escreveu: Hi, Your example isn't fully reproducible. Or maybe it just has some little errors. I can't get past this spot... ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) coordinates(ants.res) <- ~ x + y You have requested that raster::rasterToPoints() return a spatial points object, so there shouldn't be a need for sp::coordinates() to cast to a spatial object. Also, I think you are invoking "antscounts" and "antscount" as the same variable. Maybe? It's not clear to me. Finally, when you call raster::adjacent() you haven't defined the value for the cells argument when computing e1 and e2. What do you want them to be? I wish I could be more helpful. Cheers, Ben On Mon, Nov 9, 2020 at 12:46 PM ASANTOS via R-sig-Geo < r-sig-geo@r-project.org> wrote: Dear r-sig-geo Members, I'd like to find any way to create 1 (total 9 pixels) and 2 pixels (total 25 pixels) surrounding the neighborhood of each pixel (ant) in my plot image (antscount). In my example: #Packages library(spatstat) library(raster) #Selection of ants data set data(ants) geo.form<-cbind(x=ants$x,y=ants$y) #Definition of raster resolution - 10 units ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) coordinates(ants.res) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(ants.res) <- TRUE #Rasterize antscount<- rasterize(geo.form, raster(ants.res), fun='count', background=0) values(antscount)[values(antscount) > 0] = 1 #Vizualize plot(antscounts) Now, the selection of neighborhood pixels sounds easy, something like: # For 1 pixel neighborhood neigh1 <- matrix(1L, nrow=3, ncol=3) e1<-adjacent(antscounts, cells , directions=neigh1, pairs=FALSE) ng_coords1 <- xyFromCell(antscounts, e1) # For 2 pixel neighborhood neigh2 <- matrix(1L, nrow=5, ncol=5) e2<-adjacent(antscounts, cells , directions=neigh2, pairs=FALSE) ng_coords5 <- xyFromCell(antscounts, e2) But for the combination of all the information (0 and 1 new pixel values) and the new raster representation (antscounts) I don't have success. Please, any ideas? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa
[R-sig-Geo] Create pixels neighborhood in a raster
Dear r-sig-geo Members, I'd like to find any way to create 1 (total 9 pixels) and 2 pixels (total 25 pixels) surrounding the neighborhood of each pixel (ant) in my plot image (antscount). In my example: #Packages library(spatstat) library(raster) #Selection of ants data set data(ants) geo.form<-cbind(x=ants$x,y=ants$y) #Definition of raster resolution - 10 units ants.w<-as.owin(ants) ext <- as(extent(c(ants.w$xrange,ants.w$yrange)), "SpatialPolygons") ants.res<-rasterToPoints(raster(ext, resolution = 10), spatial = TRUE) coordinates(ants.res) <- ~ x + y # coerce to SpatialPixelsDataFrame gridded(ants.res) <- TRUE #Rasterize antscount<- rasterize(geo.form, raster(ants.res), fun='count', background=0) values(antscount)[values(antscount) > 0] = 1 #Vizualize plot(antscounts) Now, the selection of neighborhood pixels sounds easy, something like: # For 1 pixel neighborhood neigh1 <- matrix(1L, nrow=3, ncol=3) e1<-adjacent(antscounts, cells , directions=neigh1, pairs=FALSE) ng_coords1 <- xyFromCell(antscounts, e1) # For 2 pixel neighborhood neigh2 <- matrix(1L, nrow=5, ncol=5) e2<-adjacent(antscounts, cells , directions=neigh2, pairs=FALSE) ng_coords5 <- xyFromCell(antscounts, e2) But for the combination of all the information (0 and 1 new pixel values) and the new raster representation (antscounts) I don't have success. Please, any ideas? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] raster: extract() based in a pixel value
Dear R-sig-geo Members, I'd like to use de extract() function using the raster package for calculate a proportion of pixel with "1"s values inside a buffer given some coordinates in a raster. I try to create a function for this without success, in my hypothetical example: #Package library(raster) # Create a raster ras <- raster(ncol=1000, nrow=1000) set.seed(0) values(ras) <- runif(ncell(ras)) values(ras)[values(ras) > 0.5] = 1 values(ras)[values(ras) < 0.5] = NA # Create some coordinates pts<-sampleRandom(ras, size=30, xy=TRUE) pts.df<-as.data.frame(pts) pts.df$area<-rnorm(30, mean=10)## Here just for create a artificial covariate without any direct implication in my question #Function for extract proportion of 1s values percentual_1s<- function(x,...) { � leng1<-length(values(x) ==1) # length of 1s pixels values � lengtotal<-length(x) # total length of pixels inside buffer � perc<-(leng1/lengtotal)*100 � return(perc) } # Extract the desirable proportion in a circular 10 units buffer cent_max <- extract(ras, # raster layer ��� cbind(pts.df$x,pts.df$y),��� # SPDF with centroids for buffer ��� buffer = 10, # buffer size ��� fun=percentual_1s,�� # what to value to extract ��� df=TRUE) Here doesn't work, despite the code look like Ok. My perfect output is: #��� x� y layer� area�� percentual_1s #1 -109.26 -43.65 1 10.349010�� 23.15 #2�� 93.42 -87.21 1� 9.861920�� 45.18 #3�� 57.06� 86.85 1� 8.642071�� 74.32 #4 -109.98 -45.63 1 10.376485�� 11.56 #5� -92.34� 37.89 1 10.375138�� 56.89 #6�� 19.62� 21.51 1� 8.963949�� 88.15 Please any ideas or any another package that help in this operation? Thanks in advanced, -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- [[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] Filtering a set of points in a "ppp" object by distance using marks
Hi Marcelino, Thanks so much, I just make a little change in your code: ddd <- nndist(insects.ppp, by=factor(insects.ppp$marks)) subset(insects.ppp, marks=="termiNests" & ddd[,"antNests"] >20) I put `ddd[,"antNests"] >20` despite `ddd[,"termiNests"] >20` because I need "termiNests" mark far 20 units to each "antNests". Best wishes, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 16/09/2020 03:18, Marcelino de la Cruz Rot escreveu: 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, ___ 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 Marcelino, Thanks, I just make a little change in your code: ddd <- nndist(insects.ppp, by=factor(insects.ppp$marks)) subset(insects.ppp, marks=="termiNests" & ddd[,"antNests"] >20) I put `ddd[,"antNests"] >20` despite `ddd[,"termiNests"] >20` because I need "termiNests" mark far 20 units to each "antNests" Best wishes, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 16/09/2020 03:18, Marcelino de la Cruz Rot escreveu: 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, ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Filtering a set of points in a "ppp" object by distance using marks
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, -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- [[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] Distance rule that filtering a set of points in "ppp" class by minimum or maximum distance
Dear R-Sig_Geo Members, I'd like to find a more simple way to filtering a set of points in a "ppp" object by minimum or maximum distance. In my example using a ants ("ppp" object) in spatstat package: #Packages library(spatstat) library(sp) #Point process example data(ants) str(ants) #- attr(*, "class")= chr "ppp" #Computes the matrix of distances between all pairs of ants nests and selection of more than 20 m distance neighborhood ants.df<-as.data.frame(ants) coordinates(ants.df) <- ~x+y ants.dmat <- spDists(ants.df) #Calculate a distance matrix using spDists min.dist <- 20 #distance range ants.dmat[ants.dmat <= min.dist] <- NA # Na for less than 20 meters Here yet I need to start a relative long code for identified columns with NA in ants.dmat object, create a new data frame and than create a new ppp object based in the coordinates and marks etc... Is there any more easy way for create a new "ppp" object just only specified a distance rule that filtering a set of points in a original "ppp" object by minimum or maximum distance? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- [[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] Count occurrences less memory expensive than superimpose function in several spatial objects [SOLVED]
Thanks Vijay, Now works and I don't to need pass my shapefiles objects to ppp for using the superimpose function. The solution now is: #Packages library(spatstat) library(dplyr) library(sp) library(rgdal) library(raster) library(data.table) #Point process example data(ants) ants.df<-as.data.frame(ants) #Convert to data frame # Sample 75% in original dataset, repeat this 9 times and create a shapefile in each loop for(i in 1:9){ s.ants.df<-sample_frac(ants.df, 0.75) s.ants<-ppp(x=s.ants.df[,1],y=s.ants.df[,2],window=ants$window)#Create new ppp object sample.pts<-cbind(s.ants$x,s.ants$y) pts.sampling = SpatialPoints(sample.pts) UTMcoor.df <- SpatialPointsDataFrame(pts.sampling, data.frame(id=1:length(pts.sampling))) writeOGR(UTMcoor.df, ".",paste0('sample.shape',i), driver="ESRI Shapefile",overwrite=TRUE) } #Read all the 9 shapefiles created all_shape <- list.files(pattern="\\.shp$", full.names=TRUE) all_shape_list <- lapply(all_shape, shapefile) # target_sub1 <- rbindlist(lapply(all_shape_list, as.data.table)) res1 <- target_sub1[, .(res=.N), by=.(coords.x1,coords.x2)] res.xy1 = res1[target_sub1, on=c("coords.x1","coords.x2")] all.equal(res.xy1, res.xy1, check.attributes=FALSE) # should return TRUE #Occurrences in the same coordinate > 5 res_F<-res.xy1[res.xy1$res>5,] res_F<-as.data.frame(res_F) #Final shapefile final.pts<-cbind(res_F[,1],res_F[,2]) pts.final = SpatialPoints(final.pts) UTMcoor.df <- SpatialPointsDataFrame(pts.final, data.frame(id=1:length(pts.final))) UTMcoor.df2 <-remove.duplicates(UTMcoor.df) writeOGR(UTMcoor.df2, ".", paste0('final.ants'), driver="ESRI Shapefile", overwrite=TRUE) # Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 19/08/2020 20:49, Vijay Lulla escreveu: > Hi Alexandre, > As far as I can tell (mostly from reading the docs...no prior experience of > using multiplicity or superimpose myself) it appears that they are just > calculating the number of unique values for a combination of x,y coordinate > pairs. So, you can do this by using the group by semantics of either > tidyverse or SQL to generate the res.xy data.frame. Below is an example of > generating res.xy alternatively using data.table (I'm not as familiar with > tidyverse): > > target_sub1 <- rbindlist(lapply(target, as.data.table)) > res1 <- target_sub1[, .(res=.N), by=.(x,y)] > res.xy1 = res1[target_sub1, on=c("x","y")] > > all.equal(res.xy, res.xy1, check.attributes=FALSE) # should return TRUE > > If you're using SQL then you just join the raw table with the grouped table > and you should get the table coordinates and occurrences. And, considering > the number of coordinates you have I recommend either data.table or SQL to > generate the final output. > HTH, > Vijay. > > > On Wed, Aug 19, 2020 at 4:22 PM ASANTOS via R-sig-Geo < > r-sig-geo@r-project.org> wrote: > >> Dear r-sig-geo Members, >> >> ??? I'll like to read several shapefiles, count occurrences in the same >> coordinate and create a final shapefile with a threshold number of >> occurrences. I try to convert the shapefiles in ppp object (because I >> have some part of my data set in shapefile and another in ppp objects) >> and applied superimpose function without success. In my synthetic example : >> >> #Packages >> library(spatstat) >> library(dplyr) >> library(sp) >> library(rgdal) >> library(raster) >> >> >> #Point process example >> data(ants) >> ants.df<-as.data.frame(ants) #Convert to data frame >> >> # Sample 75% in original dataset, repeat this 9 times and create a >> shapefile in each loop >> >> for(i in 1:9){ >> s.ants.df<-sample_frac(ants.df, 0.75) >> s.ants<-ppp(x=s.ants.df[,1],y=s.ants.df[,2],window=ants$window)#Create >> new ppp object >> sample.pts<-cbind(s.ants$x,s.ants$y) >> pts.sampling = SpatialPoints(sample.pts) >> UTMcoor.df <- SpatialPointsDataFrame(pts.sampling, >> data.frame(id=1:length(pts.sampling))) >> writeOGR(UTMcoor.df, ".",paste0('sample.shape',i), driver="ESRI >> Shapefile",overwrite=TRUE) >> } >
[R-sig-Geo] Count occurrences less memory expensive than superimpose function in several spatial objects
Dear r-sig-geo Members, ??? I'll like to read several shapefiles, count occurrences in the same coordinate and create a final shapefile with a threshold number of occurrences. I try to convert the shapefiles in ppp object (because I have some part of my data set in shapefile and another in ppp objects) and applied superimpose function without success. In my synthetic example : #Packages library(spatstat) library(dplyr) library(sp) library(rgdal) library(raster) #Point process example data(ants) ants.df<-as.data.frame(ants) #Convert to data frame # Sample 75% in original dataset, repeat this 9 times and create a shapefile in each loop for(i in 1:9){ s.ants.df<-sample_frac(ants.df, 0.75) s.ants<-ppp(x=s.ants.df[,1],y=s.ants.df[,2],window=ants$window)#Create new ppp object sample.pts<-cbind(s.ants$x,s.ants$y) pts.sampling = SpatialPoints(sample.pts) UTMcoor.df <- SpatialPointsDataFrame(pts.sampling, data.frame(id=1:length(pts.sampling))) writeOGR(UTMcoor.df, ".",paste0('sample.shape',i), driver="ESRI Shapefile",overwrite=TRUE) } #Read all the 9 shapefiles created all_shape <- list.files(pattern="\\.shp$", full.names=TRUE) all_shape_list <- lapply(all_shape, shapefile) #Convert shapefile to ppp statstat target <- vector("list", length(all_shape_list)) for(i in 1:length(all_shape_list)){ target[[i]] <- ppp(x=all_shape_list[[i]]@coords[,1], y=all_shape_list[[i]]@coords[,2],window=ants$window)} #Join all ppp objects using multiplicity target_sub<-do.call(superimpose,target) res<-multiplicity(target_sub) #Occurrences in the same coordinate > 5 res.xy<-as.data.frame(target_sub$x,target_sub$y,res) res_F<-res.xy[res.xy$res>5,] #Final shapefile final.pts<-cbind(res_F[,1],res_F[,2]) pts.final = SpatialPoints(final.pts) UTMcoor.df <- SpatialPointsDataFrame(pts.final, data.frame(id=1:length(pts.final))) UTMcoor.df2 <-remove.duplicates(UTMcoor.df) writeOGR(UTMcoor.df2, ".", paste0('final.ants'), driver="ESRI Shapefile",overwrite=TRUE) This approach works very well in this synthetic example!!! But in my real data set a have the 99 shapefiles with 10^7 coordinates and when I try to use the do.call(superimpose,target) function my 32GB RAM memory crashed. Please any ideas for how I can create a new shapefile with a criteria occurrences exposed but less memory expensive than superimpose all the objects created? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- [[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] Modifications of the methods in resample() function using raster package
Thanks Frederico, But at this moment, I'll try to reduce steps (extract() - focal() - resample()) and using only resample() function to change the number of nearest cells, if is possible. First, in "bilinear" method, I'll like to change the function for work with the different number of nearest cells and second I'll expected to create a mean method to calculate the arithmetic average of the n nearest cells too. For example, in the first case for 25 nearest cells, I'll like to do something like resample (myraster, myresolution, window=matrix (nrow=5,nol=5), method="bilinear") and in the second case resample (myraster, myresolution, window=matrix (nrow=5,nol=5), fun=mean). But, obviously this commands doesn't work but expressed my general goal. Best wishes, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- Em 23/03/2020 14:23, Frederico Faleiro escreveu: > Hi Alexandre, > > you can use extract from raster (check the examples in the help page). > Depending of your data you can create a spatial polygon that cover the > number of cells you want and define the summary function or use the > buffer argument. > > Best regards, > > Frederico Faleiro > Postdoctoral Researcher in the INCT-EECBio (https://www.eecbio.ufg.br/) > Federal University of Goiás | Brazil > RG: https://www.researchgate.net/profile/Frederico_Faleiro > > > Em seg., 23 de mar. de 2020 às 10:41, ASANTOS via R-sig-Geo > mailto:r-sig-geo@r-project.org>> escreveu: > > Dear r-sig-geo members, > > ??? I?ve like to know the possibility of make some modifications in > resample() function in raster package. First, in "bilinear" method by > default is assigns a weighted average of the four nearest cells, and > I?ll like to change for five and six nearest cells, too, is possible? > Second, is it possible too to create the mean method to calculate the > arithmetic average of the four, five, and six nearest cells? > > Thanks in advanced, > Alexandre > > -- > Alexandre dos Santos > Geotechnologies and Spatial Statistics applied to Forest Entomology > Instituto Federal de Mato Grosso (IFMT) - Campus Caceres > Caixa Postal 244 (PO Box) > Avenida dos Ramires, s/n - Vila Real > Caceres - MT - CEP 78201-380 (ZIP code) > Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 > Lattes CV: http://lattes.cnpq.br/1360403201088680 > OrcID: orcid.org/-0001-8232-6722 > <http://orcid.org/-0001-8232-6722> > ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 > <http://www.researchgate.net/profile/Alexandre_Santos10> > Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ > -- > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org <mailto: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] Modifications of the methods in resample() function using raster package
Dear r-sig-geo members, ??? I?ve like to know the possibility of make some modifications in resample() function in raster package. First, in "bilinear" method by default is assigns a weighted average of the four nearest cells, and I?ll like to change for five and six nearest cells, too, is possible? Second, is it possible too to create the mean method to calculate the arithmetic average of the four, five, and six nearest cells? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Fitting an Inhomogeneous Poisson model
Dear R-Sig-Geo members, I've like to fitting an Inhomogeneous Poisson model and I have one question about better approaches and/or steps order for this. ### Simulating Inhomogeneous Poisson process library(spatstat) cI <- rpoispp(function(x, y) {200 * x^2}, lmax = 1000) cII <- rpoispp(function(x, y) {50 * y^2}, lmax = 1000) cIII <- rpoispp(function(x,y) {100 * exp(-3*x)}, lmax = 1000) ant<-superimpose(cI=cI,cII=cII,cIII=cIII) #Fitting a Poisson model with an intensity that is log-linear in the cartesian coordinates with marks m1<-ppm(ant, ~x + y + marks) # Allows different (constant) intensity for each category #Now I fitting a Inhomogeneous Poisson model with an intensity that is log-quadratic in the cartesian coordinates #with marks m2<-ppm(ant, ~polynom(x,y,2) + marks) #Comparing the models anova(m1,m2,test="Chi") #Pr(>Chi) = 0.8497 No difference and I choose the more parsimonious model m1. #But if I make the models coefficients inspection: m1#First model #� Estimate� S.E.� CI95.lo��� CI95.hi Ztest� Zval #(Intercept)� 3.0370358 0.3075106� 2.434326133� 3.6397455�� *** 9.876199 #x��� 1.3466820 0.3453355� 0.669836907� 2.0235272�� *** 3.899634 #y��� 0.6497591 0.3338675 -0.004609147� 1.3041274 1.946159 #markscII��� -1.4190842 0.2877424 -1.983048959 -0.8551194�� *** -4.931787 #markscIII�� -0.6306268 0.2154810 -1.052961828 -0.2082918��� ** -2.926601 m2 #Second model #�� Estimate� S.E.�� CI95.lo��� CI95.hi Ztest��� Zval #(Intercept)� 3.31439533 0.6095378� 2.119723� 4.5090674�� *** 5.43755529 #x��� 0.05883185 1.5736189 -3.025405� 3.1430683 0.03738634 #y��� 0.53037435 1.5536416 -2.514707� 3.5754559 0.34137497 #I(x^2)�� 1.16253158 1.2939978 -1.373658� 3.6987207 0.89840303 #I(x * y) 0.04678026 1.1589111 -2.224644� 2.3182043 0.04036570 #I(y^2)�� 0.08747001 1.2903205 -2.441512� 2.6164518 0.06778937 #markscII��� -1.41908418 0.2877424 -1.983049 -0.8551194�� *** -4.93178655 #markscIII�� -0.63062682 0.2154810 -1.052962 -0.2082918��� ** -2.92660056 I was so astonished with the no significance of I(x^2) by Ztest and my question is whats is the better workflow for improve models of this nature: 1) Try the new models with several cartesian coordinates combinations and degrees (first), 2) Try the new models with pair-wise comparison between levels (marks) (second), 3) Comparing with null model. or another approach and/or step order? -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Vila Real Caceres - MT - CEP 78201-380 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- [[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] Circular/square sample function for extract coordinates and marks in *ppp object
Dear r-sig-geo members, I've to know if there is a function for sample for extract coordinates and marks of a marked planar point pattern object? In my example: library(spatstat) ## Using a longleaf data set longleaf #Marked planar point pattern: 584 points #marks are numeric, of storage type� �double� #window: rectangle = [0, 200] x [0, 200] metres There is the spsample(x, n, "random") for example, but this function sample random points inside longleaf$window area and doesn't work for marks. I've like a function that extract longleaf points coordinates and marks inside a circular area (eg. 1 metres) or square area (eg. 1 metes h x w) around n random sample points and a data frame or better a new marked planar point patter object as output. Is there this type of function implemented in R? Thanks in advanced, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Distrito Industrial Caceres - MT - CEP 78.200-000 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Comparing distance among point pattern events
> > > # test the hypothesis that ABd is equal to ACd > > nperm <- 999 > > permout <- data.frame(ABd = rep(NA, nperm), ACd = rep(NA, nperm)) > > # create framework for a random assignment of B and C to the existing points > > BC <- superimpose(B, C) > B.len <- npoints(B) > C.len <- npoints(C) > B.sampvect <- c(rep(TRUE, B.len), rep(FALSE, C.len)) > > set.seed(2019) > for(i in seq_len(nperm)) { > B.sampvect <- sample(B.sampvect) > B.perm <- BC[B.sampvect] > C.perm <- BC[!B.sampvect] > > permout[i, ] <- c(mean(crossdist(A, B.perm)), mean(crossdist(A, C.perm))) > } > > > boxplot(permout$ABd - permout$ACd) > points(1, mean(ABd) - mean(ACd), col="red") > > table(abs(mean(ABd) - mean(ACd)) >= abs(permout$ABd - permout$ACd)) > # FALSE TRUE > # 573 426 > > sum(abs(mean(ABd) - mean(ACd)) >= abs(permout$ABd - permout$ACd)) / nperm > # 0.4264264 > > The difference between ACd and ABd is indistinguishable from that > obtained by a random resampling of B and C. > > > Sarah > > On Fri, Nov 22, 2019 at 8:26 AM ASANTOS via R-sig-Geo > wrote: >> Dear R-Sig-Geo Members, >> >> I have the hypothetical point process situation: >> >> library(spatstat) >> set.seed(2019) >> A <- rpoispp(100) ## First event >> B <- rpoispp(50) ## Second event >> C <- rpoispp(50) ## Third event >> plot(A, pch=16) >> plot(B, col="red", add=T) >> plot(C, col="blue", add=T) >> >> I've like to know an adequate spatial approach for comparing if on >> average the event B or C is more close to A. For this, I try to make: >> >> AB<-superimpose(A,B) >> ABd<-pairdist(AB) >> AC<-superimpose(A,C) >> ACd<-pairdist(A) >> mean(ABd) >> #[1] 0.5112954 >> mean(ACd) >> #[1] 0.5035042 >> >> With this naive approach, I concluded that event C is more close of A >> that B. This sounds enough for a final conclusion or more robust >> analysis is possible? >> >> Thanks in advance, >> >> Alexandre >> [[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] Comparing distance among point pattern events
Dear R-Sig-Geo Members, I have the hypothetical point process situation: library(spatstat) set.seed(2019) A <- rpoispp(100) ## First event B <- rpoispp(50) ## Second event C <- rpoispp(50) ## Third event plot(A, pch=16) plot(B, col="red", add=T) plot(C, col="blue", add=T) I've like to know an adequate spatial approach for comparing if on average the event B or C is more close to A. For this, I try to make: AB<-superimpose(A,B) ABd<-pairdist(AB) AC<-superimpose(A,C) ACd<-pairdist(A) mean(ABd) #[1] 0.5112954 mean(ACd) #[1] 0.5035042 With this naive approach, I concluded that event C is more close of A that B. This sounds enough for a final conclusion or more robust analysis is possible? Thanks in advance, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Distrito Industrial Caceres - MT - CEP 78.200-000 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Rasterize point process using specific rule
Dear members, I've like to create a raster using my event coordinates and information about the size in this coordinates as a rule for my raster creation. First, I make: library(raster) library(spatstat) #Points coordinates in UTM xp<-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) #Now I have the size of each nest in meters (marked process) area<-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 contour for the window creation W <- convexhull.xy(xp,yp) #Create ppm object syn.ppp <- ppp(x=xp,y=yp,window=W, marks=area) # Plot the point process plot(syn.ppp) #Definition of raster resolution - 1 meter square ext <- as(extent(c(W$xrange,W$yrange)), "SpatialPolygons") syn.ras.res<-rasterToPoints(raster(ext, resolution = 1), spatial = TRUE) #Now Rasterize syn.ras<- rasterize(syn.ras.res@coords, raster(syn.ras.res), *?*) I' ve like to create some kind of rule in *?*, where pixel values in my final raster was 1 when inside de area (marks=area) and zero for outside de circular area given by points coordinates (syn.ras.res@coords) using 1 meter as resolution? Any ideas, please Thanks, -- == Alexandre dos Santos Prote??o Florestal IFMT - Instituto Federal de Educa??o, Ci?ncia e Tecnologia de Mato Grosso Campus C?ceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial C?ceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley: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
Re: [R-sig-Geo] Simple Ripley's CRS test for market point patters --- addendum.
Dr. Rolf Turner, Thanks again and your explanations as a personal spatial course for me :) You are very clear, but there are some questions yet. First, these data set is not artificial, I have real coordinates and area of nests of leaf-cutting ants in my example (in my complete data set I have more than 3,000 coordinates and areas). The problem of overlap is because we represent a ants nests area as a circle (for an easy comprehension and modeling), but in the real world this areas was elliptical and some times we have an eccentricity close to 1 (in this case I have two nests close but not in overlap condition in the field). Based on our comments, I used the CSRI test, but first I try to explore more the question about intensities and size categories (s, m, and l), for this,could I test intensities by pair-wise comparison? For the last, do you have any solution for my overlap problem, if my objective is the study of my problem as a realization of a Poisson process? Please see my code below: #Packages require(spatstat) require(sp) # Create some points that represent ant nests xp<-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) #Now I have the size of each nest (marked process) area<-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 contour for the window creation W <- convexhull.xy(xp,yp) # Class of nests size - large, medium and small acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE, labels=c("s","m","l")) #Create a ppp object syn.ppp <- ppp(x=xp,y=yp,window=W, marks=acat) #Test initial hypothesis f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category. f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity for each area category. anova(f1,f2,test="Chi") #The test?? rejects the hypothesis that the intensities are the same - p =0.002015 ** *#First question, could I test intensities by pair-wise comparison?* # Small vs medium sizes acat2<-acat levels(acat2) levels(acat2)[1]<-"s_m" levels(acat2)[2]<-"s_m" levels(acat2) syn.ppp2 <- ppp(x=xp,y=yp,window=W, marks=acat2) f3 <- ppm(syn.ppp2 ~ marks) f4 <- ppm(syn.ppp2) anova(f4,f3,test="Chi") #Intensities of small and medium size are the same - p=0.237 # Medium vs large sizes acat3<-acat levels(acat3) levels(acat3)[2]<-"m_l" levels(acat3)[3]<-"m_l" levels(acat3) syn.ppp3 <- ppp(x=xp,y=yp,window=W, marks=acat3) f5 <- ppm(syn.ppp3 ~ marks) f6 <- ppm(syn.ppp3) anova(f5,f6,test="Chi") #Intensities of medium and large sizes are the different - p=7.827e-05 *** Finally small and medium sizes has the same intensity but is different of large ants size!!! With this condition I will test the CSRI: *#Now testing CSRI by simulation* foo <- function(W){ s_m <- runifpoint(length(area<55),win=W) l <- runifpoint(length(area>=55),win=W) superimpose(s_m=s_m,l=l) } simex <- expression(foo(W)) #and then set.seed(12) E <- envelope(syn.ppp2,simulate=simex,nsim=999, savefuns=TRUE) plot(E) dtst <- dclf.test(E) plot(dtst) mtst <- mad.test(E) plot(mtst) # now gives p-value = 0.003 for dtst and p-value = 0.002 for mtst Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteo Florestal IFMT - Instituto Federal de Educao, Ci??ncia e Tecnologia de Mato Grosso Campus C??ceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial C??ceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ==
Re: [R-sig-Geo] Simple Ripley's CRS test for market point patters
Thanks for your help Marcelino e for the careful explanation Rolf Turner and so sorry about my post script configuration, I expected that I solved that in my new post. First my variable area is a marked point (attribute or auxiliary information about my point process - page 7 and not a spatial covariate, effect in the outcome of my experimental area - page 50). Based in this information, my hypophyses is that the *size of ant nests* a cause of ecological intraspecific competition for resources (such as food and territory) *have different patterns of spatial distribution*, for this: #Packages require(spatstat) require(sp) # Create some points that represents ant nests xp<-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) #Now I have the size of each nest (marked process) area<-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 for the window creation W <- convexhull.xy(xp,yp) # Class of nests size - large, medium and small acat <- cut(area,breaks=c(-Inf,25,55,Inf),right=FALSE, labels=c("s","m","l")) #Create a ppp object syn.ppp <- ppp(x=xp,y=yp,window=W, marks=acat) # Test intensity hypothesis f1 <- ppm(syn.ppp) # Same (constant) intensity for each area category. f2 <- ppm(syn.ppp ~ marks) # Allows different (constant) intensity for each area category. anova(f1,f2,test="Chi") #0.002015 ** OK, the hypothesis that the intensities are the same was reject, but intensities is not the question. Based in my initial hypothesis, I've like to show envelopes and observed values of the use of some function for the three point patters (large, medium and small ant nest size)under CSR. And for this I try to: kS<-envelope(syn.ppp[syn.ppp$acat=="s"], nsim=99,fun=Kest) plot(kS,lwd=list(3,1,1,1), main="") kM<-envelope(syn.ppp[syn.ppp$acat=="m"], nsim=99,fun=Kest) plot(kM,lwd=list(3,1,1,1), main="") kL<-envelope(syn.ppp[syn.ppp$acat=="l"], nsim=99,fun=Kest) plot(kL,lwd=list(3,1,1,1), main="") But doesn't work yet. My approach now sounds correct? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteo Florestal IFMT - Instituto Federal de Educao, Ci??ncia e Tecnologia de Mato Grosso Campus C??ceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial C??ceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ == Em 23/07/2019 22:00, Rolf Turner escreveu: > > Thanks to Marcelino for chiming in on this.?? I should have responded > earlier, but was "busy on a personal matter". > > To add to what Marcelino said: > > (1) Your post in HTML was horribly garbled and a struggle to read. > *PLEASE* set your email client so that it does *NOT* post in HTML when > posting to this list. > > (2) A very minor point:?? Your construction of "syn.ppp" was > unnecessarily complicated and convoluted.?? You could simply do: > > ?? syn.ppp?? <- ppp(x=xp,y=yp,window=W,marks=area) > > (3) You appear to be confused as to the distinction between "marks" > and "covariates". These are superficially similar but conceptually > completely different.?? What you are dealing with is *marks*.?? There are > no covariates involved.?? See [1], p. 147. > > (4) Your calls to envelope() are mal-formed; the expression "area >= > 0" and "area < 25" will be taken as the values of "nrank" and ... who > knows what??? Did you not notice the warning messages you g
[R-sig-Geo] adehabitat: Working with enfa() function using *asc files
Dear R-Sig-Geo Members, ?? Many years ago I used enfa() function of adehabitat package with my *asc files with environmental variables and my script works very well, but now adehabitat is deprecated and there are many versions of adehabitat. I install the adehabitatHS but for doesn't work with my *asc files. ?? My objective was create a data2enfa object and for this I used manually all the functions in an old adehabitat package (kasc2df, import.asc, as.kasc and data2enfa) but doesn't work too. In my script I make: #Create enfa analysis object pinkdata <- data2enfa(as.kasc(list(tmax=import.asc("tmax.asc") , tmin=import.asc("tmin.asc"))), pink.sub.pnt@coords) # ??Error in count.points(pts, kasc) : ?? w should inherit the class SpatialPixels ?? I need to create a data2enfa object using *asc files as tab and organism position pink.sub.pnt@coords as pr in enfa1 <- enfa(dudi.pca(tab), pr,?? scannf = FALSE). ?? Please, any member that work with *asc files and ENFA help me? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteo Florestal IFMT - Instituto Federal de Educao, Ci??ncia e Tecnologia de Mato Grosso Campus C??ceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial C??ceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Parallel processing with extract()/randomForest() in VM
Dear R-Sig-Geo Members, ?? I create a virtual machine (VM) in Google Cloud with Ubuntu 18.04 with 8 CPU and 30 RAM memory and R 3.6.0 version, but I try to improve my spatial analysis without success or same a more faster process. If I use packages snow and doMC with all the 8 CPU's in an operation, but it use in only 12,54% of our capacity, when the objective is user extraction() in raster and RF with randomForest(). The gain of 18 seconds, I think that is not so good, then my question is there are any way for improve that? In my test, I make: # Take in the ubuntu terminal the number of processors foresteyebrazil@superforettech1:~$cat/proc/cpuinfo|grepprocess|wc-l 8 #Packages library(raster) library(snow) library(doMC) library(randomForest) registerDoMC() #Take a raster for worldclim r<-getData('worldclim', var='alt', res=5) # 1) Use extract()/ randomForest() in Virtual Machine start_time<-Sys.time() # SpatialPolygons cds1<-rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)) cds2<-rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) polys<-spPolygons(cds1, cds2) # Extract exr<-raster::extract(r, polys) tr<-ifelse(exr[[2]]<10,c("A"),c("B")) df<-cbind(tr,exr[[2]], sqrt(exr[[2]])) df2<-data.frame(as.factor(df[,1]),as.numeric(as.character(df[,2])),as.numeric(as.character(df[,3]))) df2<-df2[complete.cases(df2),] colnames(df2)<-c("res1","var1","var2") res<-NULL for(win1:9){ mod_RF<-randomForest(x=cbind(df2$var1,df2$var2), y=df2$res1, ntree=100, mtry=2) res=rbind(res,cbind(w,mean(mod_RF$err.rate[,1])*100)) } end_time<-Sys.time() end_time-start_time # #Time difference of 38.72528 secs # 2) Use extract() with snow and doMC packages in Virtual Machine start_time<-Sys.time() # SpatialPolygons cds1<-rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)) cds2<-rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) polys<-spPolygons(cds1, cds2) # Extract beginCluster(n=8) exr<-raster::extract(r, polys) tr<-ifelse(exr[[2]]<10,c("A"),c("B")) df<-cbind(tr,exr[[2]], sqrt(exr[[2]])) df2<-data.frame(as.factor(df[,1]),as.numeric(as.character(df[,2])),as.numeric(as.character(df[,3]))) df2<-df2[complete.cases(df2),] colnames(df2)<-c("res1","var1","var2") endCluster() res<-NULL mod_RF2<-foreach(1:9) %dopar%{ randomForest(x=cbind(df2$var1,df2$var2), y=df2$res1, ntree=100, mtry=2) } res=rbind(res,cbind(mean(mod_RF2$err.rate[,1])*100)) } end_time<-Sys.time() end_time-start_time # #Time difference of 20.57027 secs Thanks in advanced, -- == Alexandre dos Santos Proteo Florestal IFMT - Instituto Federal de Educao, Ci??ncia e Tecnologia de Mato Grosso Campus C??ceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial C??ceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes:http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate:www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley: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
[R-sig-Geo] Pixel calculation using extract() in geoTiff raster: faster and free of computer lock in R
Dear R-sig-geo Members, I've like to use R and raster package (extract() function) for pixel calculations (SD, Skewness and Kurtosis in 1 unit radius of a buffer) in a geoTiff raster, but my original geoTiff image has 6244(nrow), 8721(ncol) and 54453924 (ncell) dimensions. This image size causes slow computation process or computer lock, for example in my code below: ### #Packages library(raster) library(moments) #Measures of Skewness and Kurtosis ## Create artificial raster for my geoTiff simulation - dimensions of my original geoTiff: 6244, 8721, 54453924 (nrow, ncol, ncell) r <- raster(nc=8721, nr=6244) r <- setValues(r, round(runif(ncell(r))* 255)) #Extract all pixel coordinates in raster coord_r<-coordinates(r) #Extract standard deviation, skewness and kurtosis Buffer<-1 SD<-function (x, na.rm = TRUE) { if (is.matrix(x)) apply(x, 2, sd, na.rm = na.rm) else if (is.vector(x)) sqrt(var(x, na.rm = na.rm)) else if (is.data.frame(x)) sapply(x, sd, na.rm = na.rm) else sqrt(var(as.vector(x), na.rm = na.rm)) } desv_pad_R<-extract(r, coord_r, buffer = Buffer, fun = SD) str(desv_pad_R) sk_R <-extract(r,coord_r,buffer=Buffer, fun=skewness, na.rm = TRUE) str(sk_R) k_R <-extract(r,coord_r,buffer=Buffer, fun=kurtosis, na.rm = TRUE) str(k_R) # There are different approaches (eg. integration with SAGA GIS or GRASS, convert image to ASCII) for my problem? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ 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 Kernel image results in *tif format
Thank you very much Florian, Solve my problem!! Best wishes, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ == Em 27/09/2018 04:01, Florian Betz escreveu: Dear Alexandre, instead of converting the image to a SpatialPixelDataFrame converting to a raster object might be an alternative. r_pines<-raster(d_pines) writeRaster(r_pines, "Pines.tif") Regards, Florian Am 26.09.2018 um 22:25 schrieb ASANTOS via R-sig-Geo: Dear R-sig-geo Members, I've like to create Kernel image results as *tif using an object of density() function output in spatstat package. But in my example, doesn't work when I try: #Packages library(spatstat) library(raster) library(rgdal) #Swedishpines's data set in spatstat package data(swedishpines) plot(swedishpines) #CSR with K-Ripley test csr_pines <- envelope(swedishpines, Kest, nsim=99) plot(csr_pines) # r=0.75 is outside CSR #Kernel representation using 0.75 as bandwidth d_pines<-density(swedishpines, bw=0.75) plot(d_pines) #Create TIFF image r_pines <- as(d_pines, "SpatialPixelsDataFrame") writeGDAL(r_pines, "Pines.tif") # r_pines <- as(d_pines, "SpatialPixelsDataFrame") Error in as(d_pines, "SpatialPixelsDataFrame") : no method or default for coercing “im” to “SpatialPixelsDataFrame” Please any ideas for corrected this? Thanks in advanced, Alexandre ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Create Kernel image results in *tif format
Dear R-sig-geo Members, I've like to create Kernel image results as *tif using an object of density() function output in spatstat package. But in my example, doesn't work when I try: #Packages library(spatstat) library(raster) library(rgdal) #Swedishpines's data set in spatstat package data(swedishpines) plot(swedishpines) #CSR with K-Ripley test csr_pines <- envelope(swedishpines, Kest, nsim=99) plot(csr_pines) # r=0.75 is outside CSR #Kernel representation using 0.75 as bandwidth d_pines<-density(swedishpines, bw=0.75) plot(d_pines) #Create TIFF image r_pines <- as(d_pines, "SpatialPixelsDataFrame") writeGDAL(r_pines, "Pines.tif") # > r_pines <- as(d_pines, "SpatialPixelsDataFrame") Error in as(d_pines, "SpatialPixelsDataFrame") : no method or default for coercing “im” to “SpatialPixelsDataFrame” Please any ideas for corrected this? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley: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
[R-sig-Geo] Identify hotspots (centroid/geometric center when CSR distance is not satisfied) using K-Ripley approach
Dear R-sig-geo Members, I've like to identify hotspots points (centroid/geometric center of distances(r) when CSR is not satisfied), in my study case, centroids with points around 0.75 radius. This thinking in the map representation, for this objective I make: #Package library(spatstat) library(sp) library(cluster) library(lattice) #Swedishpines's data set in spatstat package data(swedishpines) plot(swedishpines) #CSR with K-Ripley test csr_pines <- envelope(swedishpines, Kest, nsim=99) plot(csr_pines) # r=0.75 is outside CSR assumption, than: ##Create matrix distance of all points coords<-cbind(swedishpines$x,swedishpines$y) res<-spDists(coords) res <- data.frame(res) # Cluster 0.75m distances clusters <- as.hclust(agnes(res, diss = T)) coords$group <- cutree(clusters, h=0.75) ## Radius 0.75 # #Visualization of centroids with points around 0.75 radius xyplot(x~y, group = group, data = coords) points(swedishpines$x,swedishpines$y, pch=16) # Doesn't work, please any ideas or new approaches? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ 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 xy coordinates from raster of interesting neighborhood cells
Thanks so much Vijay and Ben, now it works!!! Originally, I expect 24 for each 4 point (96 total). Apologies, in my example preparation I changed the raster objects 'land_mask' and 'r' and the idea of my second point [45.12069, -88.369745] is to force NA results. The code below is a solution that I expected: ### Alternative to what Ben suggested! ## Set neighborhood matrix. Focal cell is 0 while neighbors are 1!! neigh <- matrix(1L, nrow=5, ncol=5); neigh[3,3] <- 0L e1<-adjacent(r, N_cells , directions=neigh, pairs=FALSE) ng_coords1 <- xyFromCell(r, e1) plot(r); points(pts, col='red', pch=16); points(ng_coords, col='black', pch=16); points(ng_coords1, col='blue', pch=16) ### END -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ == Em 30/05/2018 16:29, Vijay Lulla escreveu: > Or you can try something like this? > > ### START > require(raster) > require(sp) > > r <- raster(nc=30, nr=30) > r <- setValues(r, round(runif(ncell(r))* 255)) > > xd <- c(-24.99270,45.12069,99.40321,73.64419) > yd <- c(-45.435267,-88.369745,-7.086949,44.174530) > pts <- data.frame(xd,yd) > pts_s<- SpatialPoints(pts) > > N_cells <- cellFromXY(r, pts_s) > e<-adjacent(r, N_cells , directions='bishop', pairs=FALSE) > > ng_coords <- xyFromCell(r, e) > plot(r) > points(pts, col='red', pch=16) > points(ng_coords, col='black', pch=16) > > [[elided Yahoo spam]] [[elided Yahoo spam]] > neigh <- matrix(1L, nrow=5, ncol=5); neigh[3,3] <- 0L > e1<-adjacent(r, N_cells , directions=neigh, pairs=FALSE) > ng_coords1 <- xyFromCell(r, e1) > plot(r); > points(pts, col='red', pch=16); > points(ng_coords, col='black', pch=16); > points(ng_coords1, col='blue', pch=16) > ### END > > HTH, > Vijay. > > On Wed, May 30, 2018 at 4:22 PM, Ben Tupper <mailto:btup...@bigelow.org>> wrote: > > Hi, > > I'm not sure why you expect 24 points - you have 4 locations and > for each you want the 4 bishop's directions - so, I think at most > you should expect 16 points. See > http://www.lpc.uottawa.ca/publications/moransi/moran.htm > <http://www.lpc.uottawa.ca/publications/moransi/moran.htm> > <http://www.lpc.uottawa.ca/publications/moransi/moran.htm > <http://www.lpc.uottawa.ca/publications/moransi/moran.htm>> for > where I expect 4 points for bishop's for each input location; > perhaps you have a different idea of bishop's direction? > > Your second point [45.12069, -88.369745] is on to the southern > edge of the raster - so it can't find points southward - only the > 2 northward. > > Finally, it may not make any difference at all, but your pts_s has > no coordinate reference. Below I show where I assign the same > projection as your raster. I also don't have 'land_mask' so I > reused 'r' and used xyFromCell to backwards extract the adjacent > cell coordinates. > > Does that do the trick for you? > > Cheers, > Ben > > ### START > require(raster) > require(sp) > > ## Create a raster > r <- raster(nc=30, nr=30) > r <- setValues(r, round(runif(ncell(r))* 255)) > plot(r) > > ##Given interesting points coordinates > xd <- c( -24.99270, 45.12069, 99.40321, 73.64419) > yd <- c(-45.435267, -88.369745, -7.086949, 44.174530) > pts <- data.frame(xd,yd) > pts_s<- SpatialPoints(pts) > projection(pts_s) <- projection(r) > points(pts_s, col="red", pch=16) > > ## Find pixels center of each point > N_cells <- cellFromXY(r, pts_s) > > e <- adjacent(r, N_cells , directions='bishop', id=TRUE, sorted = > TRUE) > xy <- xyFromCell(r, e[,'to'], spatial = TRUE) > > > #Visualization > > plot(r) > points(pts_s,
[R-sig-Geo] Extract xy coordinates from raster of interesting neighborhood cells
Dear R-Sig-Geo Members, I've like to extract xy coordinates from raster of 24 neigborhood cells that surround the 4 given points (pts_s) pixels and for this I used: #Packages require(raster) require(sp) ## Create a raster r <- raster(nc=30, nr=30) r <- setValues(r, round(runif(ncell(r))* 255)) plot(r) ##Given interesting points coordinates xd <- c(-24.99270,45.12069,99.40321,73.64419) yd <- c(-45.435267,-88.369745,-7.086949,44.174530) pts <- data.frame(xd,yd) pts_s<- SpatialPoints(pts) points(pts_s, col="red", pch=16) ## Find pixels center of each point (My idea is not to use the original point coordinates, but coordinates of the center of interesting pixel). N_cells <- cellFromXY(r, pts_s) # Extract xy coordinates from raster of 24 neighborhood cells given pixel number N_cells e<-adjacent(r, N_cells , directions='bishop', id=TRUE) coordinates(land_mask)[e[,1],] ## Doesn't return 24 neighborhood by N_cells object ng_coords<-coordinates(land_mask)[e[,1],] # #Visualization plot(r) points(pts_s, col="red", pch=16) points(ng_coords, col="black", pch=16) ## But I don't have success because the ng_coords object is not a 24 neighborhood cells of each point that I search. There are solution for this? Thanks in advance, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ 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] Estimate of the intensity of a point process, as a function of a marked point process as covariate.
Thanks a lot Dr. Turner, I will read it. Best wishes, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ == Em 31/03/2018 19:06, Rolf Turner escreveu: On 01/04/18 09:36, ASANTOS via R-sig-Geo wrote: Dear R Sig Geo Members, I've like to know if there are any function in any package for estimation density in a marked point process (e.g. geographic position and size of ants nests in square meters). My goal will be represents the density of ants nest estimated, but use nests sizes as covariate, this is possible? The Smooth.ppp() function in the spatstat package might be what you are looking for. Using nest sizes as a *covariate* does not make sense to me. A covariate should be defined at all points of the observation window, not just at the points of the observed pattern. See Adrian Baddeley, Ege Rubak, Rolf Turner (2015). Spatial Point Patterns: Methodology and Applications with R. London: Chapman and Hall/CRC Press, 2015. URL http://www.crcpress.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/9781482210200/ particularly section 5.6.4, pages 147, 148. cheers, Rolf Turner ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Estimate of the intensity of a point process, as a function of a marked point process as covariate.
Dear R Sig Geo Members, I've like to know if there are any function in any package for estimation density in a marked point process (e.g. geographic position and size of ants nests in square meters). My goal will be represents the density of ants nest estimated, but use nests sizes as covariate, this is possible? Thanks in advanced, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Create a georeferenced PDF with shapefile, points of interest and values in R
Dear Members, I created a georeferenced PDF using a shapefile and OK no problem, but I'd like to know if is possible to insert a title inside a output pdf (in top position "My geo PDF"), points of interest (pts.sampling) and values too. In my example: #Packages |require(rgdal)require(maptools)| #Create 2 polygons |sr <-SpatialPolygons(list(Polygons(list(Polygon(cbind(c(180114,180553,181127,181477,181294,181007,180409,180162,180114),c(332349,332057,332342,333250,333558,333676,332618,332413,332349,'1'),Polygons(list(Polygon(cbind(c(180042,180545,180553,180314,179955,179142,179437,179524,179979,180042),c(332373,332026,331426,330889,330683,331133,331623,332152,332357,332373,'2')))| #Convert in spatial polygon |srdf=SpatialPolygonsDataFrame(sr,data.frame(row.names=c('1','2'),PIDS=1:2))srdf@data proj4string(srdf)<-CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs")| #Create shapefile |writeOGR(srdf,getwd(),'ARS','ESRI Shapefile')| #Read shapefile |contorno_line<-readShapeLines ("ARS.shp")| #Plot |plot(contorno_line,main="My geo PDF")| #Points of interest and values |x<-c(180684.2,179786.2,180766.4,180335.5,181392.6,180881.9,181167.9,180544.5,180680.9,180259.6)y<-c(332264.8,331057.3,332190.5,331643.7,84.4,333444.8,81.6,332607.9,332625.9,331867.7)Class<-c(5,5,3,1,2,3,5,5,2,2)pts<-cbind(x,y,Class)pts.sampling =SpatialPoints(pts,proj4string=CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"))points(pts.sampling,col="red")text(pts.sampling,labels=paste(pts.sampling$Class),cex=0.7,pos=3)| #Create a geoPDF |writeOGR(contorno_line,dsn ="mygeoPDF",layer ="contorno_geo",driver ="PDF")| But my problem is, how I can put|main="My geo PDF"|, points(|pts.sampling, col="red"|) and|text(pts.sampling,labels=paste(pts.sampling$Class), cex= 0.7,pos=3)|inside the mygeoPDF.pdf created? -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes:http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 - ResearcherID: A-5790-2016 Researchgate:www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley: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
[R-sig-Geo] Extract mean pixels values using circular samples
Dear R-sig Member I'd like to extract mean pixels values using x,y coordinates by defining a buffer or area surrounding individual point coordinate. My circular sample represents the radius of a circular region around each point. I try unsuccessfully this using extract function with raster package, in my hypothetical example: require(raster) require(sp) require(plotrix)##Only for draw circular areas ## Create a raster r <- raster(nc=30, nr=30) r <- setValues(r, round(runif(ncell(r))* 255)) plot(r) #Selection of center of coordinates of my circular samples xd <- c(-24.99270,45.12069,99.40321,73.64419) yd <- c(-45.435267,-88.369745,-7.086949,44.174530) pts <- data.frame(xd,yd) pts_s<- SpatialPoints(pts) #Representation of target areas for sampling for (i in 1:length(xd)) { draw.circle(pts$xd[i],pts$yd[i],20,border="blue",col="grey") } points(pts_s, col="red", pch=16) #Extract mean pixel values in 4 circular sample with radius=20 r_sample <- extract(x = r, y = pts_s, buffer=20, fun=mean, df=TRUE) r_sample # It doesn,t work!! Could someone please help me? Thanks, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635 Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/ ___ 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] Calculate each polygon percentage inside a circles
Dear Rolf Turner, Amazing! I will never put R in doubt. My problem was solved with success, thanks very much again!!! Best wishes, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: https://www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: https://br.linkedin.com/in/alexandre-dos-santos-87961635 == Em 25/05/2016 18:09, Rolf Turner escreveu: > == > > library(spatstat) > > sr1 <- owin(poly=cbind(c(180114, 180553, 181127, 181477, 181294, > 181007, 180409, 180162, 180114), >c(332349, 332057, 332342, 333250, 333558, > 333676, 332618, 332413, 332349))) > > sr2 <- owin(poly=cbind(rev(c(180042, 180545, 180553, 180314, 179955, > 179142, 179437, 179524, 179979, 180042)), >rev(c(332373, 332026, 331426, 330889, 330683, > 331133, 331623, 332152, 332357, 332373 > > sr3 <- owin(poly=cbind(rev(c(179110, 179907, 180433, 180712, 180752, > 180329, 179875, 179668, 179572, 179269, > 178879, 178600, 178544, 179046, 179110)), >rev(c(331086, 330620, 330494, 330265, 330075, > 330233, 330336, 330004, 329783, 329665, > 329720, 329933, 330478, 331062, 331086 > > sr4 <- owin(poly=cbind(c(180304, 180403,179632,179420,180304), >c(332791, 333204, 333635, 333058, 332791))) > > wins <- solist(sr1,sr2,sr3,sr4) > > W <- union.owin(wins) > > set.seed(42) > X <- rpoispp(28/area.owin(W),win=W) > N <- npoints(X) > plot(X,cols="blue") > AD <- area.owin(disc(radius=600)) > > pct <- matrix(nrow=N,ncol=4) > rownames(pct) <- paste("point",1:N,sep=".") > colnames(pct) <- paste("sr",1:4,sep=".") > for(i in 1:npoints(X)) { > Di <- disc(radius=600,centre=c(X$x[i],X$y[i])) > for(j in 1:4) { > Aij <- intersect.owin(Di,wins[[j]],fatal=FALSE) > pct[i,j] <- if(is.null(Aij)) 0 else 100*area.owin(Aij)/AD > } > } > == [[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] [FORGED] Calculate each polygon percentage inside a circles
Dear Rolf Turner, It's much better a clean code with a minimum packages, thank you very much for your answer. But "pct" object give me a total polygon percentage around each point and I need too an identification (in columns) of individual contribution of each polygon. In my simulation, I find 50.38001% for the point 1, but this is a total percentage of polygons and I'd like to know a percentage contribution for each polygon (e.g. ID1 = 0.0 + ID1 = 10.0 + ID3 = 0.0 + ID4 = 40.38001 = 50.38001 total), this is possible? Thanks again, Best wishes, Alexandre #== library(spatstat) bdry <- list(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 332618, 332413, 332349)), cbind(rev(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 179524, 179979, 180042)), rev(c(332373, 332026, 331426, 330889, 330683, 331133, 331623, 332152, 332357, 332373))), cbind(rev(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110)), rev(c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 329783, 329665, 329720, 329933, 330478, 331062, 331086))), cbind(c(180304, 180403,179632,179420,180304), c(332791, 333204, 333635, 333058, 332791))) W <- owin(poly=bdry) set.seed(42) X <- rpoispp(28/area.owin(W),win=W) plot(X,cols="blue") for(i in 1:npoints(X)) plot(disc(radius=600,centre=c(X$x[i],X$y[i])), add=TRUE,col=rgb(1,0.5,0,alpha=0.4),border=NA) # To be fair, calculate the area of the discs using area.owin() # rather than as pi*600^2. AD <- area.owin(disc(radius=600)) pct <- numeric(npoints(X)) for(i in 1:npoints(X)) { Di <- disc(radius=600,centre=c(X$x[i],X$y[i])) pct[i] <- 100*area.owin(intersect.owin(Di,W))/AD } #== -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: https://www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: https://br.linkedin.com/in/alexandre-dos-santos-87961635 == Em 24/05/2016 23:55, Rolf Turner escreveu: > > Dear Alexandre, > > Your tasks can all be done very simply in spatstat. The code follows. > Note that I had to reverse the point order for sr2 and sr3 because > spatstat requires that the vertices of (exterior) boundaries of > polygonal windows be given in anticlockwise order. > > I commented out the plotting of the discs centred at each point of the > simulated pattern since I found the plot to be unenlightening and > messy looking. > > The resulting percentages that you are after are in "pct". > > If you want more precision for the disc areas, set npoly equal to a > large number than the default 128 (e.g.1024 or 2048) in the calls to > disc(). > > cheers, > > Rolf Turner > [[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] Calculate each polygon percentage inside a circles
Dear members, I will try to calculate each polygon percentage inside a circles given an arbitrary radius in a shapefile object with the code below and my output needs to be (Two first columns with center os circle coordinates and values of each polygon percentage): "pts$x" "pts$y" "ID1" "ID2" "ID3" "ID4" 180486.1 330358.8 16.3 0.2 NA 17.3 179884.4 331420.8 88.3 NA 96.3 NA 179799.6 333062.5 25.3 22.3 0.5 NA ... For this I try to do: #Packages require(shapefiles) require(rgdal) require(maptools) require(spatstat) require(berryFunctions) #Create 4 polygons (I create polygons because is more simple that given a shapefile in a post) sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 332618, 332413, 332349,'1') sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 331133, 331623, 332152, 332357, 332373,'2') sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110), c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 329783, 329665, 329720, 329933, 330478, 331062, 331086,'3') sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304), c(332791, 333204, 333635, 333058, 332791,'4') #Spatial Polygons sr=SpatialPolygons(list(sr1,sr2,sr3,sr4)) srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4)) #Create shapefile writeOGR(srdf, getwd(), 'POLY', 'ESRI Shapefile') #Read shapefile contorno_line_X <- readShapeSpatial("POLY.shp") plot(contorno_line_X) #Create ~28 random points in my window contorno_line_spat <- as(contorno_line_X,"owin") pts <- rpoispp(lambda=0.05, win=contorno_line_spat) points(pts, col="blue") #Set the radius for the plots radius <- 600 # radius in meters #Draw the circles using berryFunctions package for(i in 1:length(pts$x)) circle(pts$x[i],pts$y[i],r=radius, col=rgb(1,0.5,0,alpha=0.4), border=NA) # Now I'd like to calculate de polygon percentage around each point in a 600 meters radius. This is possible? There are any function for this? Thanks, Alexandre -- == Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) alexandre.san...@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 Researchgate: https://www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: https://br.linkedin.com/in/alexandre-dos-santos-87961635 == [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo