即使在转换shapefile的CRS以匹配栅格的CRS之后，也无法屏蔽空间对象

n7taea2i 于 2023-06-19 发布在其他

关注(0)|答案(2)|浏览(87)

我试图使用R中的shapefile来掩蔽光栅，特别是使用terra包。然而，我遇到了一个问题，掩码后的输出都是NaN值。我按照文档进行了操作，并确保光栅和shapefile的坐标参考系统（CRS）匹配，但问题仍然存在。
通过观察这两个图，我相信问题的原因是光栅的延伸不正确。纬度应该从-180开始到180。我使用的光栅和shapefile可以在下面的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.

来源：https://stackoverflow.com/questions/76476464/unable-to-mask-spatial-object-even-after-transforming-the-crs-of-the-shapefile-t

2条答案

按热度按时间

pqwbnv8z1#

以下是（第一个）玛格丽特的回答的简化和浓缩版本。

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)

赞(0）回复(0）举报 2023-06-19

dy2hfwbg2#

先试着旋转一下。

rotate {terra} -沿着经度旋转数据 *

将经度坐标为0到360的SpatRaster旋转到-180到180度之间的标准坐标（反之亦然）。经度在0到360之间经常用于全球气候模型。
你可能也想按shp裁剪。

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

创建于2023-06-14带有reprex v2.0.2

赞(0）回复(0）举报 2023-06-19

我来回答

即使在转换shapefile的CRS以匹配栅格的CRS之后，也无法屏蔽空间对象

2条答案

相关问题

热门标签

最新问答