1
votes

Nombre de points d'échantillonnage dans une cellule de la grille

Calculez le nombre total de points d'échantillonnage dans chaque cellule de grille d'une grille spatiale.

Je voudrais créer une grille et calculer le nombre total de points d'échantillonnage dans chaque cellule de la grille. J'ai créé des données et une grille générées aléatoirement et essayé de calculer le nombre d'enregistrements dans une grille de cellules en utilisant à la fois les packages sf et raster, en utilisant des questions SO similaires précédentes, mais sans succès. J'ai également examiné la fonction d'extraction. Je suis assez nouveau dans l'analyse spatiale.

 library(sf)
library(raster)
library(tidyverse)
library(mapview)
library(mapedit)

#Trial with sf package 
# load some spatial data. Administrative Boundary
#https://stackoverflow.com/questions/41787313/how-to-create-a-grid-of-       spatial-points
aut <- getData('GADM', country = 'aut', level = 0)
aut <- st_as_sf(aut)
#Try with polygons
grid <- aut %>% 
 st_make_grid(cellsize = 0.5, what = "polygons") %>% 
  st_intersection(aut)                               

#fake data
 lat<-runif(1000, 46.5, 48.5)
 lon<-runif(1000, 13,16)
pos<-data.frame(lat,lon)

 ggplot() + 
  geom_sf(data = aut) + 
  geom_sf(data = grid)+
geom_point(data=pos, aes(lon, lat)) 
#how to count number of records within each cell?  
 ########################################
#Trial with raster package
#https://stackoverflow.com/questions/32889531/r-how-can-i-count-how-   many-points-are-in-each-cell-of-my-grid
 r<-raster(xmn=13, ymn=46.5, xmx=16, ymx=48.5, res=0.5)
r[] <- 0
#How do I use the pos data here
xy <- spsample(as(extent(r), 'SpatialPolygons'), 100, 'random')
tab <- table(cellFromXY(r, xy))
r[as.numeric(names(tab))] <- tab
plot(r)
points(xy, pch=20)
d <- data.frame(coordinates(r), count=r[])

Je voudrais obtenir un tableau avec le nombre de points d'échantillonnage.

r sf r-raster

2 commentaires

Avez-vous examiné stat_bin_2d? ou stat_hex_bin? ceux-ci calculent le nombre et la densité par bac.

Merci pour vos conseils. Je n'avais pas envisagé stat_bin_2d ou stat_hex_bin. Comment obtenez-vous le nombre de bacs? J'ai juste commencé à jouer avec ces fonctions mais pour l'instant je n'en ai pas encore tiré quelque chose de significatif. ggplot (pos, aes (lon, lat)) + geom_bin2d (binwidth = 2) + stat_bin2d (geom = "text", aes (label = ..count ..), binwidth = 10) + scale_fill_gradient (low = "white" , haut = "rouge")

3 Réponses :

1
votes

Essayez

ggplot(pos, aes(x = lon, y=lat)) + 
geom_bin2d(binwidth = 2) +
stat_bin_2d(aes(label=stat(count)), binwidth=2, geom="text", position="identity") +
scale_fill_gradient(low = "white", high = "red")

2 commentaires

Merci. Pourriez-vous nous donner des conseils sur la façon de superposer cette figure sur une carte, par exemple la carte de l'Autriche que je donne à titre d'exemple. C'est un objet sf, mais il est possible de visialiser les objets sf avec geom_sf. J'ai essayé sans succès.

Je ne suis pas familier avec les objets sf mais chaque fois que j'ai besoin de tracer des données spartiales, je suis arrivé à la brochure ( rstudio .github.io / leaflet / forms.html )

2
votes

compter les longueurs de st_intersects (note: not st_intersection ) vous donnerait un vecteur de points contenus dans chaque cellule de la grille:

grid = st_sf(n = lengths(tab), geometry = st_cast(grid, "MULTIPOLYGON"))

mapview(grid, zcol = "n")

Si vous souhaitez ensuite lier cela à la grille en tant qu'objet sf , vous pouvez faire:

library(sf)
library(raster)
library(tidyverse)
library(mapview)
library(mapedit)

#Trial with sf package 
# load some spatial data. Administrative Boundary
#https://stackoverflow.com/questions/41787313/how-to-create-a-grid-of-       spatial-points
aut <- getData('GADM', country = 'aut', level = 0)
aut <- st_as_sf(aut)
#Try with polygons
grid <- aut %>% 
  st_make_grid(cellsize = 0.5, what = "polygons") %>% 
  st_intersection(aut)                               

#fake data
lat<-runif(1000, 46.5, 48.5)
lon<-runif(1000, 13,16)
pos<-data.frame(lat,lon)

pos = st_as_sf(pos, coords = c("lon", "lat"), crs = 4326)

tab = st_intersects(grid, pos)
lengths(tab)
[1]  0  0  0  0  4 24 23 34 23 13 14  0  0  0  0  0  0  0  3 38 40 48 46 47 33  0  0  0  0  0  0  0
[33]  0 35 48 51 35 38 44  0  0  0  0 44 43 41 53 44 32  0  0  0  0  8  8 10 12  7  0  0  0  0  0

0 commentaires

1
votes

Exemple de données

x <- rasterToPoints(r)
z <- cbind(cell=cellFromXY(r, x[,1:2]), value=x[,3])
head(z)
#     cell value
#[1,]   22    4
#[2,]   23   45
#[3,]   24   36
#[4,]   25   52
#[5,]   26   35
#[6,]   27   38

Solution

r <- rasterize(pos, r, fun="count")
plot(r)

Pour obtenir une table, vous pouvez faire

library(raster)   
aut <- getData('GADM', country = 'aut', level = 0)
r <- raster(aut, res=0.5)
lat <- runif(1000, 46.5, 48.5)
lon <- runif(1000, 13,16)
# note that you should use (lon, lat), in that order!
pos <- data.frame(lon, lat)

0 commentaires