projectRaster返回所有NA值

我正在使用來自NARCCAP的溫度.nc文件。這些數據具有極地立體投影。從最低溫度和最高溫度，我創建了一個符合楓糖漿生產日的矩陣。projectRaster返回所有NA值

我想把這個矩陣變成一個光柵，並把這個光柵投影到一個lon/lat投影。

## This is the metadata for the projection from the .nc file: 

# float lat[xc,yc] 
#   long_name: latitude 
#   standard_name: latitude 
#   units: degrees_north 
#   axis: Y 
# float lon[xc,yc] 
#   long_name: longitude 
#   standard_name: longitude 
#   units: degrees_east 
#   axis: X 
# float tasmax[xc,yc,time] 
#    coordinates: lon lat level 
#    _FillValue: 1.00000002004088e+20 
#    original_units: K 
#    long_name: Maximum Daily Surface Air Temperature 
#    missing_value: 1.00000002004088e+20 
#    original_name: T_MAX_GDS5_HTGL 
#    units: K 
#    standard_name: air_temperature 
#    cell_methods: time: maximum (interval: 24 hours) 
#    grid_mapping: polar_stereographic 

# grid_mapping_name: polar_stereographic 
# latitude_of_projection_origin: 90 
# standard_parallel: 60 
# false_easting: 4700000 
# false_northing: 8400000 
# longitude_of_central_meridian: 263 
# straight_vertical_longitude_from_pole: 263 

# The production days matrix I've created is called from a saved file: 
path.ecp2 <- paste0("E:/all_files/production/narccap/GFDL/Production_Days_SkinnerECP2", 
       year, ".RData") 
file.ecp2 <- get(load(path.ecp2)) 
dim(file.ecp2) 
# 147 116 
rast.ecp2 <- raster(file.ecp2) 
rast.ecp2 <- flip(t(rast.ecp2), 2) 
# class  : RasterLayer 
# dimensions : 116, 147, 17052 (nrow, ncol, ncell) 
# resolution : 0.006802721, 0.00862069 (x, y) 
# extent  : 0, 1, 0, 1 (xmin, xmax, ymin, ymax) 
# coord. ref. : NA 
# data source : in memory 
# names  : layer 
# values  : 0, 671 (min, max) 

# I assign the polar stereographic crs to this production days raster: 
crs("+init=epsg:3031") 
ecp2.proj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=4700000 +y_0=8400000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" 
crs(rast.ecp2) <- crs(ecp2.proj) 

rast.ecp2 
# class  : RasterLayer 
# dimensions : 116, 147, 17052 (nrow, ncol, ncell) 
# resolution : 0.006802721, 0.00862069 (x, y) 
# extent  : 0, 1, 0, 1 (xmin, xmax, ymin, ymax) 
# coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=4700000 +y_0=8400000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory 
# names  : layer 
# values  : 0, 671 (min, max)

當我使用以前爲我工作的步驟（見here），rast.ecp2的所有值去NA。我哪裏錯了？

# The projection I want to project TO: 
source_rast <- raster(nrow=222, ncol=462, xmn=-124.75, xmx=-67, ymn=25.125, ymx=52.875, 
         crs="+proj=longlat +datum=WGS84") 
rast.ecp2LL <- projectRaster(rast.ecp2, source_rast) 

rast.ecp2LL 
# class  : RasterLayer 
# dimensions : 222, 462, 102564 (nrow, ncol, ncell) 
# resolution : 0.125, 0.125 (x, y) 
# extent  : -124.75, -67, 25.125, 52.875 (xmin, xmax, ymin, ymax) 
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory 
# names  : layer 
# values  : NA, NA (min, max)

來源

2016-07-13 em.d.kear

我發佈瞭解決方案，我發現工作。它基於this post and answer。我必須首先將.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-15 17:51:12

projectRaster返回所有NA值

回答

相關問題