projectRaster gibt alle NA-Werte zurück

Question

Ich arbeite mit Temperatur .NC-Dateien von NARCCAP. Diese Daten haben eine polare stereographische Projektion. Aus Temperaturminima und -maxima habe ich eine Matrix von Tagen erstellt, die als Ahornsirupproduktionstage gelten.

Ich möchte diese Matrix in ein Raster verwandeln und dieses Raster auf eine lon/lat Projektion projizieren.

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

Wenn ich die Schritte verwenden, die vorher für mich gearbeitet (siehe here), die Werte von rast.ecp2 gehen alle auf NA. Wo gehe ich falsch?

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

Answer 1

Ich poste die Lösung, die ich gefunden habe zu arbeiten. Es basiert auf this post and answer. Ich musste zuerst die xc- und yc-Koordinaten der .nc-Datei in Längen- und Breitengrade umwandeln. Dann Ich könnte das Raster ordnungsgemäß neu projizieren. Unten ist der Code, der funktioniert hat.

Beachten Sie, dass mycrs das CRS ist, das mit der .NC-Datei "kam". Es muss dem SpatialPoints zugewiesen werden, da beim Konvertieren von xc/yc in SpatialPoints das zugehörige CRS gelöscht wird.

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) 

Von hier gehe ich auf die Matrix der Produktionstage auf die gespeicherten Dateien max.lonlat und min.lonlat basieren.

Danke! Ich hoffe, das ist hilfreich.