Как правильно настроить проекцию на растр?

Question

Я пытаюсь построить данные об осадках с метеорологического радара.Файл данных имеет матрицу 900x900 точек (900x900км).Информация о проекции из исходного файла каппи:

<projection lat_lr="48.133400" lat_ul="56.186500" type="aeqd" lon_lr="25.157600" size_x="900" size_y="900" lon_ul="11.812900">
    <lon_0>19.092600</lon_0>
    <lat_0>52.346800</lat_0>
    <ellps>+ellps=sphere</ellps>
</projection>

Я читаю файл данных (пример: https://meteomodel.pl/examples/out.txt) в матрицу и преобразую в растр:

a1 = as.matrix(read.table("/home/user/out.txt", header=F, as.is=TRUE))
a1[a1==0] <- NA
maxDBz <- 95.5
minDBz <- -31.5
step <- (maxDBz - minDBz) / 254
a1 <-  minDBz + (a1 * step)
r <- raster(a1)

ЗатемЯ пытаюсь установить экстент и CRS:

e <- extent(11.812900, 25.157600, 48.133400, 56.186500)
r <- setExtent(r, e)
crs(r) <- "+proj=aeqd +lat_0=52.346800 +lon_0=19.092600 +x_0=900 +y_0=900 +ellps=sphere +datum=WGS84 +units=km +no_defs"

Данные построены, однако проекция неверна:

https://meteomodel.pl/examples/Rplot01.png

Правильное изображение из Польского института метеорологииУправление водными ресурсами:

https://meteomodel.pl/examples/cappi.png

Что я делаю не так?

Answer 1

То, что вы делаете неправильно, это установка экстента с использованием lon / lat crs, тогда как данные имеют "+proj=aeqd. Эти должны соответствовать. Я не знаю, каков правильный размер, но вы можете приблизить его так:

p <- "+proj=aeqd +lat_0=52.346800 +lon_0=19.092600 +x_0=900 +y_0=900 +ellps=sphere +datum=WGS84 +units=km +no_defs"

e <- extent(11.812900, 25.157600, 48.133400, 56.186500)
r <- raster()
extent(r) <- e
rr <- projectExtent(r, p)
extent(rr)
#class      : Extent 
#xmin       : -541.0182 
#xmax       : 452.2122 
#ymin       : -488.8849 
#ymax       : 431.1854 

В предоставленном txt-файле указано, что желаемый экстент равен

e <- extent(-449997.470, 451000.522, -451003.637, 449998.274) 

И это говорит о том, что единицы в ваших crs должны быть m, а не km

p <- "+proj=aeqd +lat_0=52.346800 +lon_0=19.092600 +x_0=900 +y_0=900 +ellps=sphere +units=m "