Лучший способ минимизировать потерю данных при преобразовании растра из проекции в географическую c систему координат - PullRequest
Новая
       99

Лучший способ минимизировать потерю данных при преобразовании растра из проекции в географическую c систему координат

0 голосов
/ 06 августа 2020

У меня есть всемирный растровый файл в PCS примерно с разрешением 1 ° (степень составляет от 90 до -60 градусов широты и от 180 до -180 долготы), который я хочу преобразовать в GCS с разрешением 2 °, но текущий метод, который я использование приводит к искажению данных. Исходный растр содержит категориальные значения, но тот же вопрос применим к непрерывным значениям. Интересно, есть ли лучший способ сделать это. 

#create raster in PCS - CEA projection
x <- raster(ncol=360, nrow=142, xmn=-17367529, xmx=17367529, ymn=-6356742, ymx=7348382)
projection(x) <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
values(x)=runif(51120, 1, 99)
# my target raster format
y <- raster(ncol=180, nrow=75, xmn=-180, xmx=180, ymn=-60, ymx=90)
projection(y) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
values(y)=runif(13500, 1, 99)
#downscale to match my target raster resolution
x2deg=raster::aggregate(x,fact=2)
#when I try to project I get the first error
x2deg.gcs = projectRaster(x2deg, y)

Ошибка в if (maxy == miny) {: отсутствует значение, где требуется ИСТИНА / ЛОЖЬ. Вдобавок: Предупреждающее сообщение: в rgdal :: rawTransform (projfrom, projto, nrow (xy), xy [, 1], xy [,: 48 прогнозируемых точек) не конечны 

#So I slightly cut the extent and get the desired result
extent(x2deg) <- c(xmin= -17367529, xmax= 17367529, ymin= 0.999*(-6356742), ymax= .999*(7348382))
x2deg.gcs = projectRaster(x2deg, y)

Можно ли уменьшить масштаб ncol и nrows в растре x с использованием другого коэффициента? Я знаю, что мне придется сделать некоторую неравномерную интерполяцию до go с 142 строк до 75, но, например, aggregate не принимает дроби в аргументе frac. Может быть, есть другой пакет, который справляется с этим лучше, например rgdal?

1 Ответ

1 голос
/ 06 августа 2020

Это лучше работает с terra

library(terra)
x <- rast(ncol=360, nrow=142, xmin=-17367529, xmax=17367529, ymin=-6356742, ymax=7348382)
crs(x) <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
values(x) <- 1:ncell(x)

y <- rast(ncol=180, nrow=75, xmin=-180, xmax=180, ymin=-60, ymax=90)
crs(y) <- "+proj=longlat +datum=WGS84"

p <- project(x, y)
...