微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

R 光栅`projectExtent`“不能做这个转换”

如何解决R 光栅`projectExtent`“不能做这个转换”

我遇到了这个问题。从下面的大栅格中,我想剪出 100 个框,存储为 extent 对象。范围来自不同的 CRS,因此它们的单位错误。这是裁剪到法国的 CLC 土地利用栅格:

> CLCF
class      : RasterLayer 
dimensions : 11000,12000,1.32e+08  (nrow,ncol,ncell)
resolution : 100,100  (x,y)
extent     : 3100000,4300000,2100000,3200000  (xmin,xmax,ymin,ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source     : /private/var/folders/ry/d__dq69j6w31gmbnhz7ybj013sn99p/T/Rtmpd4xShx/raster/r_tmp_2021-06-11_133707_3696_56233.Grd 
names      : U2018_CLC2018_V2020_20u1 
values     : 1,128  (min,max)

enter image description here

我的边界框的 CRS 是 wgs84

> raster::extent(ms[[ix]]$inverse_mask)
class      : Extent 
xmin       : 371250 
xmax       : 381250 
ymin       : 5838250 
ymax       : 5845250 
> raster::crs(ms[[ix]]$inverse_mask)
CRS arguments:
 +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 

遵循常见的建议,这是我所做的:

empty = raster::raster()
raster::extent(empty) <- raster::extent(ms[[ix]]$inverse_mask)

> empty
class      : RasterLayer 
dimensions : 180,360,64800  (nrow,ncell)
resolution : 27.77778,38.88889  (x,y)
extent     : 371250,381250,5838250,5845250  (xmin,ymax)
crs        : +proj=longlat +datum=wgs84 +no_defs 

> raster::projectExtent(empty,crs = raster::projection(CLCF))
Error in raster::projectExtent(empty,crs = raster::projection(CLCF)) : 
  cannot do this transformation
In addition: Warning messages:
1: In showSRID(uprojargs,format = "PROJ",multiline = "NO",prefer_proj = prefer_proj) :
  discarded datum UnkNown based on GRS80 ellipsoid in Proj4 deFinition
2: In rgdal::rawTransform(projfrom,projto,nrow(xy),xy[,1],:
  663 projected point(s) not finite

在这里做错了什么?谢谢

解决方法

感谢Transform result of `st_bbox()` to other CRS

找到答案
bb = sf::st_bbox(ms[[1]]$inverse_mask)  # sf to the rescue
            bb_ll = sf::st_bbox(
                sf::st_transform( 
                    sf::st_as_sfc(bb),crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
            )
            e = raster::extent(c(bb_ll["xmin"],bb_ll["xmax"],bb_ll["ymin"],bb_ll["ymax"]))
            OL[[ix]] <- raster::crop(CLCF,raster::extent(ms[[ix]]$inverse_mask))

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。