栅格检测和/或投影crs不正确(LCC到WGS84)

Question

我有使用光栅函数brick读取的ncdf文件。我的例子/解释集中在砖块中的一层

这里提供的数据(对于tinyupload来说太大了)：https://wetransfer.com/downloads/f21973282c3f9aa4153e9b29ab26227f20190612111133/d1f6b123571e8f7f25cdc633e27c3fa220190612111133/7bd6ec

> rt <- brick("mrros/mrros_EUR-11_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_CNRM-ALADIN63_v2_mon_200601-201012.nc", 
+             varname = "mrros")[[1]] 
> rt
class      : RasterLayer 
band       : 1  (of  60  bands)
dimensions : 453, 453, 205209  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : -6.25, 5656.25, -6.25, 5656.25  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45 
source     : C:/Users/coding/Desktop/Files/mrros/mrros_EUR-11_CNRM-CERFACS-CNRM-CM5_rcp85_r1i1p1_CNRM-ALADIN63_v2_mon_200601-201012.nc 
names      : X2006.01.16 
z-value    : 2006-01-16 
zvar       : mrros 

> 
> plot(rt)
> 

​

​

光栅砖自动检测数据的crs。唯一看起来可能丢失的是数据的单位( 12.5千米间隔的数据)

> crs(rt)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45 
> crs(rt) = paste0(crs(rt), "+units=km")
> crs(rt)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45+units=km 

我想将投影更改为WGS84

> rtPr <- projectRaster(rt, crs = CRS("+init=EPSG:4326")) 
> rtPr %>% plot

​

​

这些单位显然不是WGS84，欧洲没有.12度的纵向范围(或.5度的纬度范围)！

有没有人能建议为什么projectRaster没有在适当的网格上返回数据？我假设它与输入crs的假设有关，但我不知道要将其更改为什么值，因为如果我手动执行，这些值就是我指定的值！

Answer 1

也许这会有一点帮助。

library(raster)
r <- raster(f, var="mrros")
print(r)

这显示了

#char Lambert_Conformal[]   
#  latitude_of_projection_origin: 49.5
#  standard_parallel: 49.5
#  longitude_of_central_meridian: 10.5
#  grid_mapping_name: lambert_conformal_conic

该单元存储在其他位置，栅格无法提取该单元

#dimensions
#  units: km

栅格包将创建此项目。4字符串

x <- "+proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5"

(r)gdal改变了它(我不知道为什么它添加了"lat_2")

CRS(x)
CRS arguments:
 +proj=lcc +lat_0=49.5 +lat_1=49.5 +lon_0=10.5 +ellps=WGS84 +lat_2=45 

下面的更改让它更接近了。但我只是猜测了一下。也许问题是也有轮换？？

crs(r) = "+proj=lcc +lat_1=49.5 +lon_0=-10.5 +ellps=WGS84 +lat_2=49.5 +units=km"
plot( projectRaster(r, crs="+proj=longlat +datum=WGS84") )

最好是要求提供数据源处的标准投影描述