首页 > 解决方案 > 如何使用 R 过滤掉地图形状文件之外的数据?

问题描述

我想绘制地图形状文件内的点并消除地图外的点。这是我的数据和形状文件:数据:https://www.dropbox.com/s/fsgxe8xyouky6yq/2011_50mm_map.Rda? dl=0 非洲地图形状文件:https ://www.dropbox.com/s/6vq07xiv8lk54vg/afr_g2014_2013_0.shp?dl=0

下面是我用来绘制降雨数据的脚本,但是印度洋上发生了大量降雨。我只想显示非洲大陆和马达加斯加岛内的降雨量,它们都是非洲形状文件 (afr_g2014_2013_0.shp) 的组成部分:

#Points plot inside a map
#Loading spatial analysis packages
rm(list = ls()) #Removing environment content
library(tidyverse)
library(sf)
library(maptools)
library(gplots)
library(lubridate)
library(ggspatial)
library(ggmap)
library(gstat)
library(GISTools)
library(knitr)
library(ggthemes)
library("rnaturalearth")
library("rnaturalearthdata")

#Loading the provided Africa map shape file from where it will be stored (stored in the folder below in my case)
(Africa <- st_read("C:/Users/Xolile Ncipha/Dropbox/GIS GEN data/gendata/africa/afr_g2014_2013_0.shp"))

#Load the provided Rda data file from where it will be stored 
load("2011_50mm_map.Rda")

#Plot Precipitation as spatial points
ggplot()+
  theme_minimal()+
  theme(panel.grid = element_line(color = "grey70", size = 0.1), legend.position = "right")+
  geom_point(data = Precip_2011_50mm_map_2, aes(x = Lon, y=Lat, fill = Precip), shape=21, size=3)+
  geom_sf(data = Africa, fill=NA, size=0.3, colour="grey45")+
  scale_fill_gradientn(colours = c("navyblue", "blue", "cyan", "green", "yellow", "orange", "red"))+
  #scale_fill_gradientn(colours = palett)+ # Reversed colour order
  #coord_sf(xlim = c(16, 36), ylim = c(-36, -21), expand = FALSE)+
  xlab("Longitude (degrees)") + ylab ("Latitude (degrees)")

我得到了下面的图,其中包含地图图像之外的大量数据:https ://www.dropbox.com/s/ajdjkxrbqqhln1r/Violent_Precipitation_2011_2.png?dl=0

标签: rdictionaryfilter

解决方案


以下是我将如何使用sf包解决此问题:

我无法加载您提供的 shapefile,所以我从 naturalearth 获取了一个:

africa <- rnaturalearth::ne_countries(continent = "Africa")  %>%
  sf::st_as_sf()

tmap::qtm(africa)

非洲

然后我加载您的数据集并将其转换为简单的特征(我假设它在EPSG 4326中):

load("Desktop/2011_50mm_map.Rda")
rain <- Precip_2011_50mm_map %>%
  sf::st_as_sf(coords = c("Lon", "Lat")) 

st_crs(rain) <- 4326 # this is just a guess, should be verified!

tmap::qtm(rain)

雨

最后,我使用以下方法将降雨数据剪辑/裁剪到地图上st_intersection()

rain_on_land <- st_intersection(rain, africa)

tmap::qtm(rain_on_land)

陆地上的雨


推荐阅读