r - 具有相同 CRS 但裁剪时无输出的 Shapefile 和 Raster
问题描述
问题:绘图为空,使用 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 的情节: 情节输出
# 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 的绘图:绘图输出
# Clip
Clip.step1 <- crop(SDS8_raster, extent(Terceira))
Clip <- mask(Clip.step1, Terceira)
plot(Clip)
Clip.step1 的绘图:绘图输出
Clip 的绘图:绘图输出
所有图像都显示具有相同投影的轴。我不明白我做错了什么......我是否正确定义了 CRS 和范围?
更新:
原始 shapefile 和原始 HDF 产品的全景视图: Panoply 输出
Hijmans 将两个数据集绘制在一起:绘图输出
解决方案
调试起来很棘手,因为你混合了不同的包,而你没有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\",6378206.4,294.978698213898,"
正如您所指出的(并参见此处),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)
推荐阅读
- c++ - 如何在 Qt 中的类之间发送变量
- plotly-dash - 散点图悬停文本不适用于圆形以外的标记
- apache-kafka - Kafka:隔离级别含义
- javascript - 从预签名的 AWS S3 URL 在浏览器中流式传输 WAV 文件
- node.js - 如何访问 nodejs 中的 url 查询参数?
- websphere - Linux 和 Windows 上的 WebSphere 8.5.5 之间是否有足够的差异来保证在这两个操作系统上测试我们在 WebSphere 中的应用程序?
- jquery - 即使没有指定参数名称,也可以从 URL 获取值
- google-bigquery - 如何在 Google 的 Bigquery 中获取最频繁的值
- node.js - Mac 上的 Node.js Heroku 部署 - sh: 1: nodemon: not found / npm ERR! `nodemon fileName.js`/npm 错误!在...启动脚本失败
- javascript - Vuejs-Multiselect在启用多个选项时如何仅获取选定的选项