我正在尝试使用 netcdf 以 1 度绘制来自世界海洋图集( https://www.nodc.noaa.gov/cgi-bin/OC5/woa18/woa18.pl )的海面温度数据的年度和年代际差异文件。
我已经设法创建了一张全球地图,但我无法将年度数据叠加到它上面。我的代码中是否缺少某些内容?是否可以沿 x 和 y 轴添加纬度和经度?我可以用 1955 - 64 年代数据减去 2005 - 2017 年代数据,还是需要对数据做些什么?谢谢!
remove(list=ls())
library(rgdal)
library(raster)
library(rgeos)
library(palr)
library(rworldmap)
library(reshape2)
library(raadtools)
library(raadfiles)
library(sp)
library(spbabel)
library(dplyr)
datafile <- "/Users/Desktop/woa18_decav_t00_01.nc"
path_split <- stringr::str_split(basename(datafile), ".")
img_basename <- basename(datafile)
img_startdate_string <- substr(img_basename, 2,8)
image_startdate <- as.Date(strptime(img_startdate_string, "%Y%j"))
pprj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0, 0, 0"
target_extent = c(-180, 180, -90, 90)
g4 <- rgeos::gBuffer(SpatialPoints(cbind(0, 0), proj4string = CRS(pprj)),
width = spDists(rbind(c(-180, -90), c(-180, -90)), longlat = TRUE, segments = T) * 1000,
quadsegs = 180)
data("countriesLow", package = "rworldmap")
sptable(countriesLow) <- subset(sptable(countriesLow), y_ > -84)
w <- crop(rgeos::gBuffer(spTransform(countriesLow, CRS(pprj)), width = 0), g4)
pal <- palr::sstPal(palette = TRUE)
target <- raster(g4)
res(target) <- c(25000, 25000)
sst <- projectRaster(raster(datafile), target)
sst <- mask(sst, g4)
#png("/Users/Desktop/plot1.png", width = 300, height = 300, units = "mm", res = 400)
par(mar=c(6,6,4,4))
plot(w, col="grey", breaks = pal$breaks, asp = 1, cex.axis=1.5, cex.lab=1.5, xlim = c(xmin(w), xmax(w)), ylim = c(ymin(w), ymax(w)))
mtext(text = image_startdate, side = 3, line = 0) #need this to state the decadal difference
image(chl, col = pal$cols[-1], breaks = pal$breaks, add = TRUE)
plot(g4, add=TRUE)
plot(w, col="grey", add=TRUE, breaks = pal$breaks, asp = 1, cex.axis=1.5, cex.lab=1.5, xlim = c(xmin(w), xmax(w)), ylim = c(ymin(w), ymax(w)))
#dev.off()