Search code examples
rggplot2gisrasterr-sp

ggplot gridded SpatialPolygonsDataFrame producing odd triangles and subsetting data based on point data

Using the code below I can plot the following: This code is adapted from here

enter image description here

As you can see there are few issues with the plot. I am struggling to

  1. Remove weird lines in plot
  2. Only plot cells (grids) where there are data
  3. Plot ID (see gridSpatialPolygons$values) on top of the grid cell

I realise there are a few points to this question but I hope one solution solves all.

# Load libraries
library(sp)
library(raster)
library(ggplot2)

# Projection
wgs.84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Load data
x <- c(76.82973, 76.82972, 76.82969, 76.83076, 76.83075, 76.83071, 76.83129, 76.83126, 76.83125) 
y <- c(28.26734, 28.26644, 28.26508, 28.26778, 28.26733, 28.26507, 28.26912, 28.26732, 28.26687) 
z <- c(-56.7879, -58.22462, -58.4211, -55.75333, -58.55153, -56.38619, -56.11011, -58.17415, -59.77212)

# Create data frame
dataset <- data.frame("LONGITUDE" = x, "LATITUDE" = y, "VALUES" = z)

# Create SpatialPointsDataFrame object  
datasetSP <- SpatialPointsDataFrame(coords = dataset[,c(1,2)], data = data.frame("id" = 1:nrow(dataset), "values" = dataset$VALUES), proj4string = wgs.84)

# Extent
extentDatasetSP <-extent(datasetSP)

# Make grid options

# Cell size (map units)
xCellSizeGrid <- 0.001
yCellSizeGrid <- 0.001

# Grid
grid <- GridTopology(cellcentre.offset = c(extentDatasetSP@xmin, extentDatasetSP@ymin),
                     cellsize = c(xCellSizeGrid, yCellSizeGrid), 
                     cells.dim = c(3, 7))

# Create SpatialGrid object
gridSpatial <- SpatialGrid(grid = grid, proj4string = wgs.84)


# Convert to SpatialPixels object
gridSpatialPixels <- as(gridSpatial, "SpatialPixels")

# Convert to SpatialPolygons object
gridSpatialPolygons <- as(gridSpatialPixels, "SpatialPolygons")

# Add 'id' and 'values' to every polygon
gridSpatialPolygons$id <- 1:nrow(coordinates(gridSpatialPolygons))
gridSpatialPolygons$values <- paste("Gridvalue", 1:nrow(coordinates(gridSpatialPolygons)), sep = ":")

# Get attributes from polygons 
samplePointsInPolygons2 <- datasetSP %over% gridSpatialPolygons

ggplot(gridSpatialPolygons, aes(x = long, y = lat)) + 
  geom_polygon(color = "red") +
  geom_point(data = dataset,
             aes(x = LONGITUDE,
                 y = LATITUDE))
Solution

  • When it comes to spatial objects, ggplot2 (and tidyverse in general) seems to play nicer with sf than sp. The advice below is taken from one of the help files in the associated broom package:

    Note that the sf package now defines tidy spatial objects and is the recommended approach to spatial data. sp tidiers are likely to be deprecated in the near future in favor of sf::st_as_sf(). Development of sp tidiers has halted in broom.

    Things should be fairly straightforward after conversion to sf.

    library(dplyr)

sf::st_as_sf(gridSpatialPolygons) %>%
  filter(id %in% samplePointsInPolygons2$id) %>%       # keep only grid cells with data 
  ggplot() +
  geom_sf(colour = "red") +
  geom_sf_text(aes(label = values),                    # label cells 
               nudge_y = 0.0003, colour = "grey40") +
  geom_point(data = dataset,
             aes(x = LONGITUDE,
                 y = LATITUDE))

    plot