ポイントが海のアラゴナイト飽和レベルを表すラスターからのポイント データに国の排他的経済水域を割り当てたいと考えています。ラスターは、海の多くの緯度/経度ポイントにアラゴナイト値を与える単一のレイヤーです。各緯度経度ポイントを排他的経済水域に割り当てたい。 このサイトは座標の単一のペアに対してそれを行いますが、私は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
以前、以下のコードを使用して国をラスター ポイントに割り当てましたが、これにより、国の EEZ 内にある海の多くのポイントに対して 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 と同様の機能を実行する方法はありますか? ただし、EEZ は海にありますか?
編集:再現可能な例のいくつかのデータ
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")