Преобразование экстента растрового стека из метров в координаты широты / долготы - PullRequest
0 голосов
/ 26 апреля 2018

У меня есть растровый стек в следующем формате, который я преобразовал в NetCDF, используя этот метод .Это работает, но переменные широты и долготы находятся в «метрах», который является экстентом растрового файла, и мне нужно, чтобы они были в десятичных градусах.Файл R Grid имеет значение здесь и имеет следующий формат:

class       : RasterBrick 
dimensions  : 205, 170, 34850, 12  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 160.5, 330.5, 145.5, 350.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : \TavgM_1981.grd 
names       :   Jan.1981,   Feb.1981,   Mar.1981,   Apr.1981,   May.1981,   Jun.1981,   Jul.1981,   Aug.1981,   Sep.1981,   Oct.1981,   Nov.1981,   Dec.1981 
min values  :  1.9912137,  0.8120775,  4.0446638,  4.4135274,  6.5349769,  8.6149150,  9.9991916, 11.8400562,  9.6407796,  3.5005649,  4.1205872, -0.6106244 
max values  :   9.850221,   9.121176,  10.238524,   9.858942,  11.669445,  14.260988,  15.722292,  17.235090,  15.708690,  11.598482,  11.552235,   8.981533 
time        : Jan 1981, Feb 1981, Mar 1981, Apr 1981, May 1981, Jun 1981, Jul 1981, Aug 1981, Sep 1981, Oct 1981, Nov 1981, Dec 1981 

В идеале я хотел бы преобразовать его в файл NetCDF, который имеет следующий формат для чтения моделью:

Dimensions:    (lat: 408, lon: 881, nb2: 2, time: 660)
Coordinates:
  * lon        (lon) float64 -44.81 -44.69 -44.56 -44.44 -44.31 -44.19 ...
  * lat        (lat) float64 21.81 21.94 22.06 22.19 22.31 22.44 22.56 22.69 ...
  * time       (time) datetime64[ns] 1951-01-16T12:00:00 1951-02-16T12:00:00 ...
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) datetime64[ns] 1951-01-16T12:00:00 ...
    tas        (time, lat, lon) float64 nan nan nan nan nan nan nan nan nan ...

Я пытался преобразовать экстент Растрового кирпича в широты / долготы, используя этот метод:

sputm <- SpatialPoints(test@extent@xmin, proj4string=CRS("+proj=utm +zone=29N +datum=WGS84"))

Но это выдает эту ошибку:

Error in (function (classes, fdef, mtable)  : unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric"’

РЕДАКТИРОВАТЬ: также попытался

crs(test) <- '+proj=merc +datum=WGS84'
x <- projectRaster(test, crs='+proj=longlat +datum=WGS84') 

writeRaster(x, "rstack2.nc", overwrite=TRUE, format="CDF",     varname="Temperature", varunit="degC", 
        longname="Temperature -- raster stack to netCDF", xname="Longitude",   yname="Latitude", zname="Time (Month)")

Но когда я проверяю это в целом, координаты широты и долготы равны 0,0.

Я не уверен, куда идти дальше.

1 Ответ

0 голосов
/ 27 апреля 2018

Как только я узнал Ирландию, я посмотрел на прогнозы Ирландии.Я думаю, что это TM75 / Irish Grid.Я изменил разрешение вашего растра так, чтобы оно составляло 2500 м.

library(raster)

r <- stack("~/Bureau/TavgM_1981")
xmax(r) <- xmin(r) + 2500 * ncol(r)
ymax(r) <- ymin(r) + 2500 * nrow(r)

crs(r) <- "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +ellps=mod_airy +towgs84=482.5,-130.6,564.6,-1.042,-0.214,-0.631,8.15 +units=m +no_defs"

# Convert to degrees
x_wgs84 <- projectRaster(r, crs='+proj=longlat +datum=WGS84') 

writeRaster(x, "~/Bureau/TavgM_wgs84.tif", overwrite = TRUE)

Ниже вы увидите результат этого разрешения + изменение проекции по сравнению с целью.
Мы не так далеки от реальности, нов исходном растре есть что-то плохое.Если это (к сожалению) обычно не иметь crs, когда кто-то дает вам слой, это обычно не пропускает правильное разрешение.
Поскольку предоставленный вами растр сохраняется как файл ".gri / .grd", этоозначает, что он был создан с R. Вы должны попытаться получить оригинальный код R, который создал файл, вместо окончательного измененного слоя.

Irish raster after resolution and crs modified

...