正確にラスターイメージにマップを追加するには

Question

私は海面温度データをプロットしようとしており、データがNAsと混同されないようにカラーの土地イメージを追加しようとしています。私はそうするために複数の方法を試しましたが、下の画像に示すように、マップはデータに対して適切に並んでいません。

この問題が再現可能にするために、ここで私が働いているファイルとのDropboxへのリンクです：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、および

#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) 

をプロット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変換およびベース私はこの空間データのために使用してい

-The投影が「EPSG:4326

-Hereは私も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引数ですが、擬似メルカトル投影をオプションとして持たないようです。

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))) 

あなたがArcMapのか、QGISで適切にそれを得ることができ、あなたの目的に十分である分解能でRasterLayerを設定し、

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)