Search code examples
rgeospatialkmlr-raster

Counting number of points on a raster layer in R

I've got a map with certain number of points on it. I want to (1) calculate the number of points that fall within the raster layer, and (2) extract these points to a data frame.

This is what I've done:

# 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)

I want to calculate and extract black points located on the colorful parts of this map enter image description here

Solution

  • First raster::projectRaster does not "set" the projection but, rather reprojects the raster given a transformation and resampling. Given the computational requirements of this it is much faster to reproject the point data using sp::spTransform. Once your data is in the same projection space, you can use raster::extract to extract the raster values. Values out side the raster or in nodata (NA) areas will be assigned NA values. You can drop these observations using a simple NA index with which.

    It looks like your data may have a constant value outside of the permafrost. Once you identify what this value is (eg., 0) you can remove these points as well. Here is a worked example. First we add packages and create some example data that is similar to yours.

    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)

    Now, we extract the raster values and look at the dimensions of the point object (please note that I do not have to use the sp=TRUE argument in the raster::extract function). 

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

    Now we create a row index indicating which rows contain zeros, make sure that we have identified rows (using an if statement) and then remove them. Looking at the object dimensions again, we can see how many points were removed from the original point data. 

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