Растр обнаруживает и / или неправильно проецирует crs (LCC - WGS84) - PullRequest
Новая
       43

Растр обнаруживает и / или неправильно проецирует crs (LCC - WGS84)

0 голосов
/ 12 июня 2019

У меня есть ncdf файлы, которые я читаю, используя растровую функцию brick. Мой пример / объяснение фокусируется на одном слое в кирпиче

Доступные здесь данные (слишком велики для крошечной загрузки): https://wetransfer.com/downloads/f21973282c3f9aa4153e9b29ab26227f20190612111133/d1f6b123571e8f7f25cdc633e27c3fa220190612111133/7bd6ec

> rt <- brick("mrros/mrros_EUR-11_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_CNRM-ALADIN63_v2_mon_200601-201012.nc", 
+             varname = "mrros")[[1]] 
> rt
class      : RasterLayer 
band       : 1  (of  60  bands)
dimensions : 453, 453, 205209  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : -6.25, 5656.25, -6.25, 5656.25  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45 
source     : C:/Users/coding/Desktop/Files/mrros/mrros_EUR-11_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_CNRM-ALADIN63_v2_mon_200601-201012.nc 
names      : X2006.01.16 
z-value    : 2006-01-16 
zvar       : mrros 

> 
> plot(rt)
>

enter image description here

Растровый кирпич автоматически обнаруживает Crs данных. Единственное, что может показаться отсутствующим, - это единицы данных (данные с интервалами 12,5 км) 

> crs(rt)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45 
> crs(rt) = paste0(crs(rt), "+units=km")
> crs(rt)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45+units=km

Я хочу изменить проекцию на WGS84 

> rtPr <- projectRaster(rt, crs = CRS("+init=EPSG:4326")) 
> rtPr %>% plot

enter image description here

Эти единицы явно не WGS84, в Европе нет продольного диапазона .12 градусов (или широтного диапазона .5 градусов)!

Кто-нибудь может посоветовать, почему projectRaster не возвращает данные по соответствующей сетке? Я предполагаю, что это как-то связано с предположением о вводе crs, но я понятия не имею, на что его изменить, так как это те значения, которые я бы указал, если бы сделал это вручную!

1 Ответ

0 голосов
/ 12 июня 2019

Возможно, это немного помогает 

library(raster)
r <- raster(f, var="mrros")
print(r)

Это показывает 

#char Lambert_Conformal[]   
#  latitude_of_projection_origin: 49.5
#  standard_parallel: 49.5
#  longitude_of_central_meridian: 10.5
#  grid_mapping_name: lambert_conformal_conic

Устройство хранится в другом месте, и оно не принимается растром 

#dimensions
#  units: km

Растровый пакет создает эту строку proj.4 

x <- "+proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5"

(r) gdal меняет ее (я не знаю, почему она добавляет "lat_2") 

CRS(x)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45

Приведенные ниже изменения получаютсяэто близко.Но я просто догадываюсь немного.Возможно, проблема в том, что есть также вращение ??

crs(r) = "+proj=lcc +lat_1=49.5 +lon_0=-10.5 +ellps=WGS84 +lat_2=49.5 +units=km"
plot( projectRaster(r, crs="+proj=longlat +datum=WGS84") )

Возможно, лучше спросить стандартное описание проекции у источника данных

...