Почему обрезка растрового стека меняет названия слоев? - PullRequest
Новая
       20

Почему обрезка растрового стека меняет названия слоев?

1 голос
/ 27 февраля 2020

Ежегодно я обрабатываю многослойные файлы netCDF с данными о суточных осадках от CHIRPS. У меня есть файлы для всего мира, каждый файл размером около 1,2 ГБ. Мне нужно рассчитать индексы на основе данных об осадках для каждой ячейки в растре для определенного региона c. Чтобы сделать это, я пытаюсь обрезать файлы, чтобы получить прямоугольную форму angular над моей областью интересов, используя пакет растровых изображений.

Это код, который я использую, пример для первый файл. 

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(79, 89, 25, 31), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# Set directory with original files
setwd("~/data")

# Read file
chirps81 <- stack("chirps-v2.0.1981.days_p05.nc")
chirps81crop <-crop(chirps1981, crop_extent)

# Write cropped file back to different folder
setwd("~/croppeddata")
writeRaster(chirps81crop, "chirps81crop.nc", overwrite=TRUE)

По какой-то причине, однако, при записи файла слои теряют свое имя. В исходных файлах и после обрезки имена имеют названия слоев в формате «X1981.01.01». Но после записи и чтения файла netCDF с new file <- stack("chirps81crop.nc") имена слоев меняются на формат «X1» до «X365». Я думаю, что с этим должно быть все в порядке, если предположить, что порядок слоев не перепутался, но я не понимаю, что происходит с именами слоев, и если это происходит, потому что с кодом что-то не так.

Ответы [ 2 ]

2 голосов
/ 27 февраля 2020

Это функция writeRaster(), которая теряет имена слоев, а не операцию обрезки. Можно использовать функции более низкого уровня ncdf, чтобы назначить значение цифры c (к сожалению, не строку) каждому слою, который затем будет отображаться в названии слоев после чтения. Вдохновленный примером здесь , я создал некоторый код, который показывает это. 

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(5.74, 5.75, 50.96, 50.97), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# make a sample file
r <- raster(system.file("external/test.grd", package="raster"))
r.latlon <- projectRaster(r, crs = proj4string(crop_extent))
writeRaster(x=r.latlon, filename = 'test.nc', format = 'CDF', overwrite=TRUE)

# read the sample as a 2 layer stack and crop it
test <- stack('test.nc', 'test.nc')
writeRaster(test, 'teststack.nc', overwrite=TRUE, format='CDF')
testcrop <- crop(test, crop_extent)
names(testcrop)
# [1] "test.1" "test.2"

# write the cropped file and make the zname equal to Layer
writeRaster(testcrop, 'testcrop.nc', overwrite=TRUE, format='CDF', zname='Layer')
# open the cdf file directly
nc <- nc_open('testcrop.nc', write = T)
# give the layers numbers starting from 10 so 
# we can see them easily
layers = 1:nlayers(testcrop) + 10
layers
# [1] 11 12
ncvar_put(nc, 'Layer', layers)
nc_close(nc)

newtestcrop <- stack('testcrop.nc')
names(newtestcrop)
# [1] "X11" "X12"
nc <- nc_open('testcrop.nc', write = F)
layers = ncvar_get(nc, 'Layer')
layers
# [1] 11 12
nc_close(nc)

Так что при написании растра можно получить имена с номерами под вашим контролем, но я не не знаю достаточно о вашей среде, чтобы определить, поможет ли это, так как может быть сложно сопоставить нужные имена с одним однозначным числом.

1 голос
/ 27 февраля 2020

Надеюсь, вы не возражаете, если я предложу решение, не относящееся к R, но эта задача намного проще из командной строки с использованием CDO: 

cdo sellonlatbox,79,89,25,31 chirps-v2.0.1981.days_p05.nc cropped_file.nc

Какие показатели вы хотите рассчитать? Я подозреваю, что их можно быстро и легко рассчитать с помощью функций CDO ...

...