1

我想将国家专属经济区分配给栅格中的点数据，其中的点代表海洋中的文石饱和度。栅格是一个单一图层，可为海洋中的许多纬度/经度点提供 Aragonite 值。我想将每个纬度/经度点分配给一个专属经济区。 该站点针对单对坐标执行此操作，但我有 15,000 个点，因此我希望可以在 R 中执行此操作。

数据如下所示：

      long      lat Aragonite
1 20.89833 84.66917  1.542071
2 22.69496 84.66917  1.538187
3 24.49159 84.66917  1.537830
4 26.28822 84.66917  1.534834
5 28.08485 84.66917  1.534595
6 29.88148 84.66917  1.532505

以前我使用下面的代码将国家分配给栅格点，但这会为国家专属经济区内的海洋中的许多点返回 NA。

#convert the raster to points for assigning countries
r.pts <- rasterToPoints(r, spatial = TRUE)

#view new proj 4 string of spatialpointsdataframe
proj4string(r.pts)
##[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

###converting reclassified points to countries
# The single argument to this function, points, is a data.frame in which:
#   - column 1 contains the longitude in degrees
#   - column 2 contains the latitude in degrees
coords2country = function(r.pts)
{
countriesSP <- getMap(resolution='high')
#countriesSP <- getMap(resolution='high') #you could use high res map from rworldxtra if you were concerned about detail

#setting CRS directly to that from rworldmap
r.pts = SpatialPoints(r.pts, proj4string=CRS(proj4string(countriesSP))) 

# use 'over' to get indices of the Polygons object containing each point 
indices = over(r.pts, countriesSP)
# return the ADMIN names of each country
indices$ADMIN 
#indices$ISO3 # returns the ISO3 code
#indices$continent   # returns the continent (6 continent model)
#indices$REGION   # returns the continent (7 continent model)
}

#get country names for each pair of points
rCountries <- coords2country(r.pts)

除了海洋中的专属经济区，有没有办法对 coords2countries 执行类似的功能？

编辑：可重现示例的一些数据

dput(head(r.pts))
structure(list(layer = c(5, 5, 5, 5, 5, 5), x = c(-178.311660375408,-176.511660375408, -174.711660375408, -172.911660375408, -171.111660375408,-169.311660375408), y = c(73.1088933113454, 73.1088933113454,73.1088933113454, 73.1088933113454, 73.1088933113454, 73.1088933113454),.Names = c("layer", "x", "y"),row.names = c(NA, 6L), class = "data.frame")
4

1 回答 1

0

您需要一个包含专属经济区的 shapefile。下载此处的版本：World EEZ v9 (2016-10-21, 123 MB) http://www.marineregions.org/downloads.php#marbound

readOGR()您可以使用包中的函数加载 EEZ shapefile rgdal。将 EEZ shapefile zip 解压缩到您的工作目录并运行countriesSP <- rgdal::readOGR(dsn = ".", layer = "eez_boundaries")代替countriesSP <- getMap(resolution='high')

您提供的示例数据都不属于一个国家的专属经济区，所以我无法判断这是否真的有效，但它可能应该......

library(sp)
library(rworldmap)
library(rgeos)

r <- read.table(header = TRUE, text = "
long lat Aragonite
1 20.89833 84.66917  1.542071
2 22.69496 84.66917  1.538187
3 24.49159 84.66917  1.537830
4 26.28822 84.66917  1.534834
5 28.08485 84.66917  1.534595
6 29.88148 84.66917  1.532505
")
# or
#r <- data.frame(long = c(-178.311660375408,-176.511660375408, -174.711660375408, -172.911660375408, -171.111660375408,-169.311660375408), 
#                lat = c(73.1088933113454, 73.1088933113454,73.1088933113454, 73.1088933113454, 73.1088933113454, 73.1088933113454))

r.pts <- sp::SpatialPoints(r)

# download file from here: http://www.marineregions.org/download_file.php?fn=v9_20161021
# put the zip file in your working directory: getwd()
unzip('World_EEZ_v9_20161021.zip')

# countriesSP <- rworldmap::getMap(resolution = "high")
# or
countriesSP <- rgdal::readOGR(dsn = ".", layer = "eez_boundaries")

r.pts <- sp::SpatialPoints(r.pts, proj4string = sp::CRS(proj4string(countriesSP))) 
indices <- over(r.pts, countriesSP)
indices$ADMIN
于 2017-02-10T01:10:13.433 回答