정확하게 래스터 이미지에지도를 추가하는 방법

Question

데이터를 NAs과 혼동하지 않도록 해수면 온도 데이터를 플롯하고 토지의 컬러 이미지를 추가하려고합니다. 그렇게하기 위해 여러 가지 방법을 시도했지만 아래 이미지에서 볼 수 있듯이 맵은 데이터와 관련하여 올바르게 정렬되지 않습니다.

는 여기에 내가 함께 일하고 있어요 파일과 보관에 대한 링크입니다,이 문제를 재현 할 수 있도록하려면 다음

https://www.dropbox.com/s/e8pwgmnhvw4s0nf/sst.nc4?dl=0 내가 지금까지 개발 한 코드가;

library(ncdf4) 
library(raster) 
library(mapdata) 
library(mapproj) 
library(rgeos) 
library(ggplot2) 

비아 ncdf4, rasterToPoints, map_data 및 ggplot2

eight = nc_open("Downloads/sst.nc4") 
sst = ncvar_get(eight, "sst") 
sst = raster(sst) 
sst = t(flip(sst, 1)) # have to orient the data properly 

# extract the dimensions and set the extent 
lat.min = min(eight$dim$lat$vals) 
lat.max = max(eight$dim$lat$vals) 
lon.min = min(eight$dim$lon$vals) 
lon.max = max(eight$dim$lon$vals) 
sst = setExtent(sst, ext = c(lon.min, lon.max, lat.min, lat.max)) 

# provide proper projection 
crs(sst) = "+init=epsg:4326" 

# convert raster to points 
sst.p <- rasterToPoints(sst) 
df <- data.frame(sst.p) 
colnames(df) <- c("Longitude", "Latitude", "sst") 
usa = map_data("usa") 
ggplot(data=df, aes(y=Latitude, x=Longitude)) + 
    geom_raster(aes(fill=sst)) + 
    theme_bw() + 
    coord_equal() + 
    scale_fill_gradient("SST (Celsius)", limits=c(0,35)) + 
    geom_polygon(data = usa, aes(x=long, y = lat, group = group)) + 

    theme(axis.title.x = element_text(size=16), 
     axis.title.y = element_text(size=16, angle=90), 
     axis.text.x = element_text(size=14), 
     axis.text.y = element_text(size=14), 
     panel.grid.major = element_blank(), 
     panel.grid.minor = element_blank(), 
     legend.position = "right", 
     legend.key = element_blank() 
) 

비아 래스터가 maptools 데이터, SP 변형 및베이스

#read in the data 
sst = raster("Downloads/sst.nc4", varname = "sst", stopIfNotEqualSpaced=FALSE) 

# get world map data 
data("wrld_simpl", package="maptools") 

## Crop to the desired extent, then plot 
newext <- c(lon.min, lon.max, lat.min, lat.max) 
out <- crop(wrld_simpl, newext) 

#transform to proper CRS 
out = spTransform(out, "+init=epsg:4326") 

#plot 
plot(out, col="khaki", bg="azure2") 
plot(sst, add = T) 

0,123,516 플로팅 나는이 공간 데이터를 사용하고

년 - 투사가 'EPSG:4326

- 여기가 나는 또한 mapproj으로 map() 기능을 사용하려고 한 sst.nc4 출력 투사

<crs>PROJCS["Mercator_1SP/World Geodetic System 1984", 
     GEOGCS["World Geodetic System 1984", 
     DATUM["World Geodetic System 1984", 
     SPHEROID["WGS 84", 6378135.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
     AUTHORITY["EPSG","6326"]], 
     PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
     UNIT["degree", 0.017453292519943295], 
     AXIS["Geodetic longitude", EAST], 
     AXIS["Geodetic latitude", NORTH]], 
     PROJECTION["Mercator_1SP"], 
     PARAMETER["latitude_of_origin", 0.0], 
     PARAMETER["central_meridian", 0.0], 
     PARAMETER["scale_factor", 1.0], 
     PARAMETER["false_easting", 0.0], 
     PARAMETER["false_northing", 0.0], 
     UNIT["m", 1.0], 
     AXIS["Easting", EAST], 
     AXIS["Northing", NORTH]]</crs> 

를 지시하는 XML 조각입니다 s projection 인수이지만 옵션으로 의사 - mercator 투영이없는 것 같습니다.

Answer 1

이 부분은 다소 혼란 스럽습니다. 그래서

"cells are not equally spaced; you should extract values as points" 

의 그렇게하자 : 일반적으로 가장 쉬운 방법은이 파일에 대한, 즉,이 오류를 제공,

sst = raster("sst.nc4", varname = "sst") 

수 있지만 것

library(ncdf4) 
library(raster) 
library(maptools) 

d <- nc_open("sst.nc4") 
sst <- ncvar_get(d, "sst") 
lon <- ncvar_get(d, "lon") 
lat <- ncvar_get(d, "lat") 
nc_close(d) 

xy <- cbind(rep(lon, length(lat)), rep(lat, each=length(lon))) 

결합하고 (약 NA 값을 제거 세포의 절반 ...

xyv <- na.omit(cbind(xy, as.vector(sst))) 

는 점

r <- raster(extent(range(lon), range(lat)), res=1/6) 
r <- rasterize(xyv[, 1:2], r, xyv[,3], fun=mean) 

플롯

data(wrld_simpl) 
w <- crop(wrld_simpl, r) 

plot(r) 
plot(w, col='gray', add=TRUE)