Shapefile 和 Raster 具有相同的 CRS 但裁剪时没有输出

问题描述

问题：绘图为空，使用 shapefile 裁剪栅格时没有结果。

我使用的是带有原始 epsg4326 投影的 shapefile 和带有原始正弦投影的 MOdis 产品。正如您在脚本中看到的那样，我将两者都转换为相同的投影 (DesiredCRS)，但是在制作栅格剪辑时，我没有得到任何结果。

library(terra)
library(sf) 
library(sp) 
library(raster)

# Inputs
HDFfile <- "ModisProductsOriginal/MCD18A1.A2001001.h15v05.061.2020097222704.hdf"
Shapefile <- "Shapefile/Outline_5/Portugal_Outline_5_CAOP2019.shp"
DesiredCRS <- "+proj=longlat +datum=wgs84"

# Feature
Shp <- st_read(Shapefile)
Terceira <- Shp[Shp$DI == '43',1]
Terceira <- st_transform(Terceira,DesiredCRS)
plot(Terceira,axes=TRUE)

Terceira 的情节： Plot Output

# Modis data
SDSs <- sds(HDFfile)
SDS8 <- SDSs[8]
SDS8_template <- rast(ncol=1200,nrow=1200,xmin=-34,xmax=-24,ymin=36,ymax=41,crs=DesiredCRS)
SDS8_reprojected <- project(SDS8,SDS8_template) # Reproject changes pixels?
SDS8_raster <- raster(SDS8_reprojected)
plot(SDS8_raster)

SDS8_raster 的绘图：Plot Output

# Clip
Clip.step1 <- crop(SDS8_raster,extent(Terceira))
Clip <- mask(Clip.step1,Terceira)
plot(Clip)

Clip.step1 的情节：Plot Output

剪辑的情节：Plot Output

所有图像都显示具有相同投影的轴。我不明白我做错了什么......我是否正确定义了 CRS 和范围？

更新：

原始 shapefile 和原始 HDF 产品的全景视图： Panoply Output

Hijmans 将两个数据集绘制在一起：Plot Output

解决方法

这很难调试，因为你混合了不同的包，而且你没有show(object)。无论如何，这是一种 terra 方法，显示在何处查看对象元数据会很有用；以及同时显示栅格和矢量数据的图。

library(terra)
HDFfile <- "MCD18A1.A2001001.h15v05.061.2020097222704.hdf"
Shapefile <- "Portugal_Outline_5_CAOP2019.shp"
DesiredCRS <- "+proj=longlat +datum=WGS84"

Shp <- vect(Shapefile)
Terceira <- Shp[Shp$DI == '43',1]

SDSs <- sds(HDFfile)
SDS8 <- SDSs[8]

首先将矢量数据投影到原始栅格数据。

TercSin <- project(Terceira,crs(SDS8))
plot(SDS8,1)
lines(TercSin)

如您所见，这行不通。原因是 GDAL/PROJ 从文件中读取或假设了错误的椭球

cat(substr(crs(SDS8),1,230),"\n")
#PROJCRS["unnamed",#    BASEGEOGCRS["Unknown datum based upon the Clarke 1866 ellipsoid",#        DATUM["Not specified (based on Clarke 1866 spheroid)",#            ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,

或者像这样

describe("HDF4_EOS:EOS_GRID:\"MCD18A1.A2001001.h15v05.061.2020097222704.hdf\":MODISRAD:GMT_1200_DSR")[1:8]

#[1] "Driver: HDF4Image/HDF4 Dataset"                                         
#[2] "Files: MCD18A1.A2001001.h15v05.061.2020097222704.hdf"                   
#[3] "Size is 1200,1200"                                                     
#[4] "Coordinate System is:"                                                  
#[5] "PROJCRS[\"unnamed\","                                                   
#[6] "    BASEGEOGCRS[\"Unknown datum based upon the Clarke 1866 ellipsoid\","
#[7] "        DATUM[\"Not specified (based on Clarke 1866 spheroid)\","       
#[8] "            ELLIPSOID[\"Clarke 1866\","      

正如您指出的（并参见 here），MODIS（或至少某些产品）使用

modcrs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m"

那我们试试吧

TercSin <- project(Terceira,modcrs)
plot(SDS8,1)
lines(TercSin)

看起来没问题。所以我们需要做

crs(SDS8) <- modcrs

并继续

Terceira <- project(Terceira,DesiredCRS)
Terceira

SDS8_template <- rast(ncol=1200,nrow=1200,xmin=-34,xmax=-24,ymin=36,ymax=41,crs=DesiredCRS)
SDS8_reprojected <- project(SDS8,SDS8_template) 
SDS8_reprojected

# Again plot the two datasets together
plot(SDS8_reprojected)
lines(Terceira)

Clip.step1 <- crop(SDS8_reprojected,Terceira)
Clip <- mask(Clip.step1,Terceira)
Clip

plot(Clip)
lines(Terceira)