将雷达数据重新投影到不同的坐标系

Question

如果我使用此代码将一个 spat 栅格重新投影到另一个栅格：

old_crs<-"+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752 +units=m"
a<-terra::rast(ncol = 700, nrow = 765, ext=c(-236275.403,  106985.856, 501792.127,  900792.861), crs=old_crs)

newcrs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
b <- terra::rast(ncols=700, nrows=765, ext=c(3.2, 7.3, 50.7, 53.7), crs=newcrs)
w <- terra::project(a,b)

重新投影光栅 a:

class       : SpatRaster 
dimensions  : 765, 700, 1  (nrow, ncol, nlyr)
resolution  : 490.3732, 521.5696  (x, y)
extent      : -236275.4, 106985.9, 501792.1, 900792.9  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs 
>

到光栅 b

class       : SpatRaster 
dimensions  : 765, 700, 1  (nrow, ncol, nlyr)
resolution  : 0.005857143, 0.003921569  (x, y)
extent      : 3.2, 7.3, 50.7, 53.7  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 

并创建光栅 w

class       : SpatRaster 
dimensions  : 765, 700, 1  (nrow, ncol, nlyr)
resolution  : 0.005857143, 0.003921569  (x, y)
extent      : 3.2, 7.3, 50.7, 53.7  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 

然后它就可以正常工作了。

一旦我用这个光栅替换 Ratser a：

class       : SpatRaster 
dimensions  : 765, 700, 1  (nrow, ncol, nlyr)
resolution  : 490.3732, 521.5696  (x, y)
extent      : -236275.4, 106985.9, 501792.1, 900792.9  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs 
source(s)   : memory
name        : lyr.1 
min value   :     0 
max value   :   133 

这是 similat ro Ratser a 但有值，投影不起作用。

Error: [project] Cannot do this transformation
In addition: Warning messages:
1: In x@ptr$warp(y@ptr, "", method, mask[1], align[1], FALSE, opt) :
  GDAL Error 1: PROJ: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body

我做错了什么？

我想重新投影这个光栅

class       : SpatRaster 
dimensions  : 765, 700, 1  (nrow, ncol, nlyr)
resolution  : 490.3732, 521.5696  (x, y)
extent      : -236275.4, 106985.9, 501792.1, 900792.9  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs 
source(s)   : memory
name        : lyr.1 
min value   :     0 
max value   :   133 

到“纬度/经度”坐标。

谢谢！

Answer 1

您的数据

library(terra)
a <- terra::rast(ncol = 700, nrow = 765, ext=c(-236275.403,  106985.856, 501792.127,  900792.861), crs=old_crs)
values(a) <- 1:ncell(a)
old_crs<-"+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752 +units=m"
newcrs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

old_crs

不好。椭球参数的值以公里为单位，而不是以米为单位，这使得地球非常小。我想这就是你收到“其他天体警告”的原因。

此外，您不应该将 +ellps 和 +towgs84 添加到

new_crs

（但这不是这里的问题）。

这样比较合理：

old_crs<-"+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378137 +b=6356752 +units=m"
newcrs <- "+proj=longlat +datum=WGS84"

调试此问题的最佳方法是仅提供 crs，而不是模板栅格，以查看数据的最终位置。使用固定的 CRS，投影“有效”，但数据不在正确的位置。

crs(a) <- old_crs
w <- terra::project(a, newcrs)
plot(w)

下面使用的CRS获取的是荷兰所在区域的数据

crs(a) <- "+proj=stere +lat_0=46 +lon_0=6"
w <- terra::project(a, newcrs)
round(ext(w), 3)
#SpatExtent : 2.413, 7.624, 50.461, 54.083 (xmin, xmax, ymin, ymax)

但这不太可能是正确的。所以我认为你需要询问数据提供者并询问他们 CRS 实际上是什么。或者看看 KNMI 提供的其他一些数据集——这些数据集可能有正确的 CRS。