问问题
984 次
1 回答
1
我正在发布我发现可行的解决方案。它基于这篇文章和答案。我必须首先将 .nc 文件的 xc 和 yc 坐标转换为经度和纬度点。然后我可以正确地重新投影光栅。下面是有效的代码。
请注意,这mycrs
是 .nc 文件“附带”的 CRS。它必须分配给SpatialPoints
从 xc/yc 转换到SpatialPoints
删除相关 CRS 的原因。
years <- seq(from=1971, to=2000, by=5)
model <- "CRCM"
convert.lonlat <- function(model, year)
{
max.stem <- "E:/all_files/www.earthsystemgrid.org/CCSM/tasmax_"
inputfile <- paste0(max.stem, model, "_ccsm_", year, "010106.nc")
lat <- raster(inputfile, varname="lat")
lon <- raster(inputfile, varname = "lon")
plat <- rasterToPoints(lat)
plon <- rasterToPoints(lon)
lonlat <- cbind(plon[,3], plat[,3])
lonlat <- SpatialPoints(lonlat, proj4string = crs(base.proj))
mycrs <- crs("+proj=stere +lon_0=263 +x_0=3475000 +y_0=7475000 +lat_0=90 +ellps=WGS84")
plonlat <- spTransform(lonlat, CRSobj = mycrs)
maxs <- brick(inputfile, varname="tasmax")
projection(maxs) <- mycrs
extent(maxs) <- extent(plonlat)
max.lonlat <- projectRaster(maxs, base.proj)
save(max.lonlat, file=paste0("E:/all_files/production/narccap/CCSM/", model, "max_lonlat_", year, ".RData"))
min.stem <- "E:/all_files/www.earthsystemgrid.org/CCSM/tasmin_"
inputfile <- paste0(min.stem, model, "_ccsm_", year, "010106.nc")
lat <- raster(inputfile, varname="lat")
lon <- raster(inputfile, varname = "lon")
plat <- rasterToPoints(lat)
plon <- rasterToPoints(lon)
lonlat <- cbind(plon[,3], plat[,3])
lonlat <- SpatialPoints(lonlat, proj4string = crs(maurer.proj))
mycrs <- crs("+proj=stere +lon_0=263 +x_0=3475000 +y_0=7475000 +lat_0=90 +ellps=WGS84")
plonlat <- spTransform(lonlat, CRSobj = mycrs)
mins <- brick(inputfile, varname="tasmin")
projection(mins) <- mycrs
extent(mins) <- extent(plonlat)
min.lonlat <- projectRaster(mins, maurer.proj)
save(min.lonlat, file=paste0("E:/all_files/production/narccap/CCSM/", model, "min_lonlat_", year, ".RData"))
}
lapply(years, convert.lonlat, model=model)
从这里我继续根据保存的文件max.lonlat
和min.lonlat
.
谢谢!希望这会有所帮助。
于 2016-08-15T17:51:12.580 回答