英文:
Unable to mask spatial object even after transforming the CRS of the shapefile to match the CRS of the raster
问题
我正在尝试使用R中的terra
包来使用形状文件对栅格进行蒙版处理。但是,我遇到了一个问题,即在蒙版处理后的输出中全部是NaN
值。我已经按照文档的说明确保了栅格和形状文件的坐标参考系统(CRS)匹配,但问题仍然存在。
通过观察这两个图,我认为问题的原因是栅格的扩展不正确。纬度应该从-180开始,一直到180。我使用的栅格和形状文件可以在以下链接中获取:链接。
# ---- NOAA全球温度数据 ----
# 加载terra包
if (!require(terra)) install.packages('terra', repos='https://rspatial.r-universe.dev')
# 加载NOAA全球温度数据
tmax_2018 <- rast("data/temperature/tmax.2018.nc")
tmax_2018_Jan1 <- tmax_2018[[1]]
> tmax_2018_Jan1
class : SpatRaster
dimensions : 360, 720, 1 (nrow, ncol, nlyr)
resolution : 0.5, 0.5 (x, y)
extent : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : tmax.2018.nc
varname : tmax (每日最高温度)
name : tmax_1
unit : °C
time : 2018-01-01 UTC
# ---- 厄瓜多尔形状文件 ----
# 加载sf包
if (!require("sf")) install.packages("sf")
# 加载厄瓜多尔形状文件
ecuador_shp <- st_read('data/ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)
# 更改CRS以匹配SpatRaster的CRS
ecuador_shp <- st_transform(ecuador_shp, crs(tmax_2018_Jan1))
# 从sf类转换为terra的SpatVector类
ecuador_shp <- vect(ecuador_shp)
> ecuador_shp
class : SpatVector
geometry : polygons
dimensions : 24, 9 (geometries, attributes)
extent : -92.00854, -75.20007, -5.017001, 1.681093 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
names : ID_0 ISO NAME_0 ID_1 NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
type : <num> <chr> <chr> <num> <chr> <chr> <chr> <chr> <chr>
values : 68 ECU Ecuador 1 Azuay Provincia Province NA NA
68 ECU Ecuador 2 Bolivar Provincia Province NA NA
68 ECU Ecuador 3 Cañar Provincia Province NA NA
# ---- 蒙版处理 ----
tmax_2018_Jan1_ecu <- terra::mask(tmax_2018_Jan1, ecuador_shp)
plot(tmax_2018_Jan1) # 显示一个空白图。
英文:
I'm trying to mask a raster using a shapefile in R, specifically using the terra
package. However, I'm encountering an issue where the output after masking is all NaN
values. I've followed the documentation and made sure that the coordinate reference systems (CRS) of both the raster and shapefile match, but the problem persists.
By looking at the two plots, I believe the cause of the problem is that the extension of the raster is incorrect. The latitude should start at -180 up to 180. The raster and shapefiles I am using can be obtain in the following link.
# ---- NOAA global temperature data ----
# Load the terra package
if (!require(terra)) install.packages('terra', repos='https://rspatial.r-universe.dev')
# Load NOAA global temperature data
tmax_2018 <- rast("data/temperature/tmax.2018.nc")
tmax_2018_Jan1 <- tmax_2018[[1]]
> tmax_2018_Jan1
class : SpatRaster
dimensions : 360, 720, 1 (nrow, ncol, nlyr)
resolution : 0.5, 0.5 (x, y)
extent : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : tmax.2018.nc
varname : tmax (Daily Maximum Temperature)
name : tmax_1
unit : degC
time : 2018-01-01 UTC
# ---- Ecuador shapefile ----
# Load sf package
if (!require("sf")) install.packages("sf")
# Load Ecuador shapefiles
ecuador_shp <- st_read('data/ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)
# Change the CRS to match the SpatRaster's CRS
ecuador_shp <- st_transform(ecuador_shp, crs(tmax_2018_Jan1))
# Transforming sf class to SpatVector class from terra
ecuador_shp <- vect(ecuador_shp)
> ecuador_shp
class : SpatVector
geometry : polygons
dimensions : 24, 9 (geometries, attributes)
extent : -92.00854, -75.20007, -5.017001, 1.681093 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
names : ID_0 ISO NAME_0 ID_1 NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
type : <num> <chr> <chr> <num> <chr> <chr> <chr> <chr> <chr>
values : 68 ECU Ecuador 1 Azuay Provincia Province NA NA
68 ECU Ecuador 2 Bolivar Provincia Province NA NA
68 ECU Ecuador 3 Cañar Provincia Province NA NA
# ---- Masking ----
tmax_2018_Jan1_ecu <- terra::mask(tmax_2018_Jan1, ecuador_shp)
plot(tmax_2018_Jan1) # Shows an empty plot.
答案1
得分: 2
尝试先旋转它。
> rotate {terra} - 绕经度旋转数据
> 将具有从0到360的经度坐标的SpatRaster旋转为标准坐标,范围在-180到180度之间(或反之亦然)。在全球气候模型中经常使用0到360之间的经度。
你可能也想按shp文件裁剪。
library(sf)
#> 链接到 GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1;sf_use_s2() 为 TRUE
library(terra)
#> terra 1.7.23
# 载入,检查范围
tmax_2018 <- rast("tmax.2018.nc")
ext(tmax_2018)
#> 空间范围 : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)
# 旋转,检查范围
tmax_2018_r <- rotate(tmax_2018)
ext(tmax_2018_r)
#> 空间范围 : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
par(mfrow = c(2, 1))
plot(tmax_2018[[1]])
plot(tmax_2018_r[[1]])
<!-- -->
ecuador_shp <- read_sf('ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)
tmax_2018_crop <- crop(tmax_2018_r, ecuador_shp)
plot(tmax_2018_crop[[1]])
plot(ecuador_shp$geometry, add = TRUE)
tmax_2018_masked <- mask(tmax_2018_crop, ecuador_shp)
plot(tmax_2018_masked[[1]])
plot(ecuador_shp$geometry, add = TRUE)
<!-- -->
tmax_2018_masked
#> 类型 : SpatRaster
#> 尺寸 : 13, 34, 365 (nrow, ncol, nlyr)
#> 分辨率 : 0.5, 0.5 (x, y)
#> 空间范围 : -92, -75, -5, 1.5 (xmin, xmax, ymin, ymax)
#> 坐标参考 : 经度/纬度 WGS 84
#> 来源 : 内存
#> 名称 : tmax_1, tmax_2, tmax_3, tmax_4, tmax_5, tmax_6, ...
#> 最小值 : 9.768158, 17.85175, 15.38072, 17.51550, 13.29921, 13.41261, ...
#> 最大值 : 33.109238, 34.98264, 35.73123, 33.20931, 33.43045, 33.31937, ...
#> 时间 : 2018-01-01 到 2018-12-31 UTC
<sup>创建于 2023-06-14,使用 reprex v2.0.2</sup>
英文:
Try rotating it first.
> rotate {terra} - Rotate data along longitude
> Rotate a SpatRaster that has longitude coordinates from 0 to 360, to standard coordinates between -180 and 180 degrees (or vice-versa). Longitude between 0 and 360 is frequently used in global climate models.
And you may want to crop by shp as well.
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.7.23
# load, check extent
tmax_2018 <- rast("tmax.2018.nc")
ext(tmax_2018)
#> SpatExtent : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)
# rotate, check extent
tmax_2018_r <- rotate(tmax_2018)
ext(tmax_2018_r)
#> SpatExtent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
par(mfrow = c(2, 1))
plot(tmax_2018[[1]])
plot(tmax_2018_r[[1]])
<!-- -->
ecuador_shp <- read_sf('ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)
tmax_2018_crop <- crop(tmax_2018_r, ecuador_shp)
plot(tmax_2018_crop[[1]])
plot(ecuador_shp$geometry, add = TRUE)
tmax_2018_masked <- mask(tmax_2018_crop, ecuador_shp)
plot(tmax_2018_masked[[1]])
plot(ecuador_shp$geometry, add = TRUE)
<!-- -->
tmax_2018_masked
#> class : SpatRaster
#> dimensions : 13, 34, 365 (nrow, ncol, nlyr)
#> resolution : 0.5, 0.5 (x, y)
#> extent : -92, -75, -5, 1.5 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> names : tmax_1, tmax_2, tmax_3, tmax_4, tmax_5, tmax_6, ...
#> min values : 9.768158, 17.85175, 15.38072, 17.51550, 13.29921, 13.41261, ...
#> max values : 33.109238, 34.98264, 35.73123, 33.20931, 33.43045, 33.31937, ...
#> time : 2018-01-01 to 2018-12-31 UTC
<sup>Created on 2023-06-14 with reprex v2.0.2</sup>
答案2
得分: 2
这是margusl答案的简化和精简版本。
library(terra)
tmax_2018 <- rast("tmax.2018.nc")
tmax_2018 <- rotate(tmax_2018)
ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
tmax_2018_ecu <- crop(tmax_2018, ecu, mask=TRUE)
plot(tmax_2018_ecu, 1)
lines(ecu)
如果您不想旋转栅格数据,可以改变矢量数据的位置。
library(terra)
tmax_2018 <- rast("tmax.2018.nc")
ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
secu <- shift(secu, 360)
tmax_2018_secu <- crop(tmax_2018, secu, mask=TRUE)
plot(tmax_2018_secu, 1)
lines(secu)
然后您可以将结果移动。如果这是一个问题,这可能更有效。
tmax <- shift(tmax_2018_secu, -360)
英文:
Here is (first) a simplified and condensed version of margusl's answer.
library(terra)
tmax_2018 <- rast("tmax.2018.nc")
tmax_2018 <- rotate(tmax_2018)
ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
tmax_2018_ecu <- crop(tmax_2018, ecu, mask=TRUE)
plot(tmax_2018_ecu, 1)
lines(ecu)
If you do not want to rotate the raster data, you could shift the vector data instead.
library(terra)
tmax_2018 <- rast("tmax.2018.nc")
ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
secu <- shift(secu, 360)
tmax_2018_secu <- crop(tmax_2018, secu, mask=TRUE)
plot(tmax_2018_secu, 1)
lines(secu)
And you could then shift the result. That could be more efficient if that is an issue.
tmax <- shift(tmax_2018_secu, -360)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论