r - 如何使用 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
解决方案
以下是我将如何使用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)