1
4

1 回答 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.lonlatmin.lonlat.

谢谢!希望这会有所帮助。

于 2016-08-15T17:51:12.580 回答