R: ошибка при перемещении растра с помощью функции Shift: он НЕ двигается. Почему? - PullRequest
1 голос
/ 23 октября 2019

У меня есть файл nc (который я могу прочитать как растр в R и ArcGIS) и шейп-файл (полигон) для Канады.

Когда я открываю оба файла в ArcMap, он перепроектирует их на лету, и многоугольник будет расположен в правильном месте:

enter image description here

Затем я использую R, чтобы открыть и построить их, используя код:

## set the working directory:
setwd("C:/Users/Desktop/Data")

## required libraries
library(raster)
library(rgdal)
library(ncdf4)

## reading the polygon:
My_study_area <- readOGR(dsn = ".", layer = "Canada_AB")

## create a connection to the downloaded NC file:
fn <- "tasmin_day_CanESM5_ssp126_r1i1p1f1_gn_21010101-23001231_cliped_AB.nc"
ncin <- nc_open(fn) 
ncin

## create a raster from the nc file:
fn_brick <- brick(fn) 
fn_brick
extent(fn_brick)

Вот вывод:

class      : Extent
xmin       : 229.2188 
xmax       : 271.4062
ymin       : 39.06842
ymax       : 64.18258


## since it has 73,000 layers, let us only work on one layer: fn_brick[[1]]

Моя цель:

## I want to shift this raster in the x direction for -177.1876 units:  
## why -177.1876 units ?
## because I know if I shift it for -177.1876 units, the raster and the polygon will be in the right place 

fn_brick_shifted <- shift(fn_brick[[1]], dx= -177.1876, dy= 0, filename="new_nc_file", format="GTiff", datatype='FLT4S', overwrite=T)

extent(fn_brick_shifted)

Несмотря на смещение растра, его экстент все еще остается неизменным:

class      : Extent
xmin       : 229.2188 
xmax       : 271.4062
ymin       : 39.06842
ymax       : 64.18258

Однако, давайте построим их:

plot(fn_brick_shifted)
lines(My_study_area)

Результат вR выглядит следующим образом:

enter image description here

Как видно, растр НЕ смещен , поэтому полигон не нанесен на графикправильное место.

Однако мы можем видеть информацию о координатах для файлов:

Для растра:

crs(fn_brick_shifted)
CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

Для многоугольника:

crs(My_study_area)
CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

Я допустил какую-либо ошибку при использовании функции shift для растра?

Нужно ли использовать другую функцию для получения моей цели?

Любые комментарии или помощь будутвысоко ценится

Если вы считаете, что наличие растровых и полигональных файлов может помочь, вы можете загрузить их здесь в виде zip-файла:

nc и polygon.zip

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...