projectRaster возвращает все значения NA

Я работаю с файлами температуры .nc из NARCCAP. Эти данные имеют полярную стереографическую проекцию. Из минимальных и максимальных температур я создал матрицу дней, которые можно квалифицировать как дни производства кленового сиропа.

Я хочу превратить эту матрицу в растр и спроецировать этот растр на проекцию долгота/широта.

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

Когда я использую шаги, которые работали для меня ранее (см. ">здесь), все значения 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)

coordinates na raster projection extent
person em.d.kear    schedule 13.07.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я отправляю решение, которое я нашел для работы. Он основан на этом сообщении и ответе. Сначала мне пришлось преобразовать координаты xc и yc файла .nc в точки долготы и широты. Тогда я мог правильно перепроецировать растр. Ниже приведен код, который сработал.

Обратите внимание, что mycrs — это CRS, который «пришел с» файлом .nc. Он должен быть назначен 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.

Спасибо! Надеюсь, это полезно.

person em.d.kear    schedule 15.08.2016