I am working with temperature .nc files from NARCCAP. These data have a polar stereographic projection. From temperature minimums and maximums, I've created a matrix of days that qualify as maple syrup production days.
I want to turn this matrix into a raster, and project this raster to a lon/lat projection.
## 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)
When I use the steps that worked for me previously (see here), the values of rast.ecp2 all go to NA. Where am I going wrong?
# 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)
I am posting the solution I found to work. It is based on this post and answer. I had to first convert the xc and yc coordinates of the .nc file to longitude and latitude points. Then I could properly reproject the raster. Below is the code that worked.
Note that mycrs is the CRS that "came with" the .nc file. It has to be assigned to the SpatialPoints since converting from xc/yc to SpatialPoints drops the associated 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)
From here I go on to make the matrix of production days based on the saved files max.lonlat and min.lonlat.