2016-12-12 18 views
0

국가에서 배타적 경제 구역을 지정하여 점이 바다의 아라 고 네이트 채도 수준을 나타내는 래스터의 데이터를 가리 키려고합니다. 래스터는 바다의 많은 위도/경도 점에 대해 Aragonite 값을 제공하는 단일 레이어입니다. 배타적 경제 수역에 각 위도/경도 포인트를 지정하고 싶습니다. This site 좌표의 한 쌍을 위해 그것을 않습니다하지만 난 15,000 포인트가 나는 R.에서 할 수 있습니다 바라고 있도록국가 경계가 아닌 독점 경제 구역에서 R의 coords2country 함수 사용

데이터는 다음과 같이 :

 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") 
+1

재현성 –

+0

에 대한 샘플 데이터의'dput을()'주십시오 죄송합니다, 나는 데이터의 부분 집합으로 너무 이상 했어요. –

답변

0

에 대한 일부 데이터 당신은 EEZs을 포함하는 Shape 파일이 필요합니다. 여기에 하나를 다운로드 : 세계 EEZ의 V9 (2016년 10월 21일 123 메가 바이트) http://www.marineregions.org/downloads.php#marbound

당신은 rgdal 패키지에서 readOGR() 기능으로 EEZ의 Shape 파일을로드 할 수 있습니다. EEZ shapefile zip을 작업 디렉토리에 압축을 풀고 countriesSP <- getMap(resolution='high')

제공된 예제 데이터는 국가의 EEZ에 해당하지 않으므로 실제로 작동하는지는 알 수 없습니다. .

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