Подсчет количества точек на растровом слое в R - PullRequest
Новая
       81

Подсчет количества точек на растровом слое в R

0 голосов
/ 04 марта 2020

У меня есть карта с определенным количеством точек на ней. Я хочу (1) вычислить количество точек, попадающих в растровый слой, и (2) извлечь эти точки во фрейм данных.

Вот что я сделал: 

# Packages
library(raster)
library(ggplot2)
library(maptools)
library(tidyverse)
library(dplyr)
library(sp)


# Transform tree ring kml to dataframe
zz<-getKMLcoordinates('treering.kml', ignoreAltitude=TRUE)
l<-as.data.frame(zz)
l<-t(l)

tree <-SpatialPointsDataFrame(l, l,
                                  proj4string = CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84 
                                  +no_defs +towgs84=0,0,0"))


# Get world map
data(wrld_simpl)

# Transform World to raster
r <- raster(wrld_simpl, res = 1)
wrld_r <- rasterize(wrld_simpl, r)

# Import permafrost layer to raster
dist1<-raster("PZI.flt")


# Set CRS
dist1 <- projectRaster(from = dist1, crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
+towgs84=0,0,0"))

# Change colours

micolor <- rev(rainbow(12, alpha = 0.35))
transp <- rainbow(12, alpha = 0)
micolor[1:3] <- transp[1]


# Plot all
plot(wrld_r, col = "lightgrey")
plot(dist1, add=TRUE, legend = F, col = micolor)
plot(tree, add=T, pch = 20, col='black', cex=0.2)

Я хочу вычислить и извлечь черные точки, расположенные на красочных частях этой карты enter image description here

1 Ответ

1 голос
/ 04 марта 2020

First raster::projectRaster не "устанавливает" проекцию, а скорее перепроектирует растр с учетом преобразования и повторной выборки. Принимая во внимание вычислительные требования этого, гораздо быстрее перепроецировать данные точек, используя sp::spTransform. Как только ваши данные окажутся в том же пространстве проекции, вы можете использовать raster::extract для извлечения растровых значений. Значениям вне растра или в областях узлов (NA) будут присвоены значения NA. Вы можете отбросить эти наблюдения, используя простой индекс NA с which.

Похоже, ваши данные могут иметь постоянное значение вне вечной мерзлоты. Как только вы определите, что это за значение (например, 0), вы также можете удалить эти точки. Вот рабочий пример. Сначала мы добавляем пакеты и создаем некоторые примеры данных, которые похожи на ваши. 

library(sp)
library(raster)

dist1 <- raster(nrow=20, ncol=20)
  dist1[] <- sample(1:10, ncell(dist1), replace=TRUE)
    dist1[200:400] <- 0

trees <- sampleRandom(dist1, 100, sp=TRUE)

plot(dist1)
  plot(trees,pch=20,col="red",add=TRUE)

Теперь мы извлекаем растровые значения и просматриваем размеры точечного объекта (обратите внимание, что мне не нужно использовать аргумент sp = TRUE в функции raster::extract). 

trees@data <- data.frame(trees@data, dist1 = extract(dist1, trees))
  dim(trees)

Теперь мы создадим индекс строки, указывающий, какие строки содержат нули, убедитесь, что мы идентифицировали строки (используя оператор if), а затем удалим их. Снова посмотрев на размеры объекта, мы увидим, сколько точек было удалено из исходных данных точек. 

( idx <- which(trees$dist1 %in% 0) )
  if(length(idx) > 0) trees <- trees[-idx,]
    dim(trees)
...