英文:
Reprojecting radar data to a different coordinate system
问题
如果我使用以下代码将一个空间栅格重新投影到另一个：
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 从原始坐标系重新投影到 b 的坐标系。
如果要将类似于 a 的带有值的栅格重新投影到“Lat/Lon”坐标系，并出现错误，请检查以下内容：
- 栅格
a是否具有正确的坐标引用（+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）。 - 确保栅格
a和栅格b具有相同的分辨率和大小。 - 检查栅格
a是否有合适的投影参数，以确保与新坐标系兼容。
如果仍然出现错误，请提供更多关于栅格 a 的信息，以便更详细地解决问题。
英文:
If I use this code to reproject one spat raster to another:
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)
Which reprojects raster a:
class : SpatRasterdimensions : 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>
To Raster b
class : SpatRasterdimensions : 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
And creates Raster w
class : SpatRasterdimensions : 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
Then it works without a problem.
As soon as I replace Ratser a with this raster:
class : SpatRasterdimensions : 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_defssource(s) : memoryname : lyr.1min value : 0max value : 133
Which is similat ro Ratser a but has values, The projection does not work.
Error: [project] Cannot do this transformationIn 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
What am I doing wrong?
I would like to reporject this Raster
class : SpatRasterdimensions : 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_defssource(s) : memoryname : lyr.1min value : 0max value : 133
To "Lat/Lon" coordinates.
Thank you!
答案1
得分: 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<-"+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(a) <- old_crsw <- terra::project(a, newcrs)plot(w)
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)
请注意，以上只是代码部分的翻译。
英文:
Your data
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 is not good. The values for the ellipsoid parameters are in km, not in m, and that makes the earth very small. I suppose that is why you get the "other celestial body" warning.
Also, you should not add +ellps and +towgs84 to new_crs (but that is not the problem here).
This is more reasonable:
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"
The best way to debug this is to only provide a crs, not a template raster, to see where the data end up. With the fixed CRS the projection "works", but the data are not in the right place.
crs(a) <- old_crsw <- terra::project(a, newcrs)plot(w)
The CRS used below gets the data in the area where The Netherlands is.
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)
But this is not likely to be correct. So I think you need to ask the data provider what the CRS actually is. Or perhaps look at some other datasets that the KNMI provides --- these might have the correct CRS.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论