无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。

huangapple go评论73阅读模式
英文:

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) # 显示一个空白图。

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。

英文:

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(&#39;terra&#39;, repos=&#39;https://rspatial.r-universe.dev&#39;)

# Load NOAA global temperature data
tmax_2018 &lt;- rast(&quot;data/temperature/tmax.2018.nc&quot;)

tmax_2018_Jan1 &lt;- tmax_2018[[1]]
&gt; 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

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。

# ---- Ecuador shapefile ----
# Load sf package
if (!require(&quot;sf&quot;)) install.packages(&quot;sf&quot;)

# Load Ecuador shapefiles
ecuador_shp &lt;- st_read(&#39;data/ecu_shapefiles/ECU_adm1.shp&#39;)
ecuador_shp &lt;- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)

# Change the CRS to match the SpatRaster&#39;s CRS
ecuador_shp &lt;- st_transform(ecuador_shp, crs(tmax_2018_Jan1))

# Transforming sf class to SpatVector class from terra
ecuador_shp &lt;- vect(ecuador_shp)
&gt; 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        : &lt;num&gt; &lt;chr&gt;   &lt;chr&gt; &lt;num&gt;   &lt;chr&gt;     &lt;chr&gt;     &lt;chr&gt;     &lt;chr&gt;     &lt;chr&gt;
 values      :    68   ECU Ecuador     1   Azuay Provincia  Province        NA        NA
                  68   ECU Ecuador     2 Bolivar Provincia  Province        NA        NA
                  68   ECU Ecuador     3   Ca&#241;ar Provincia  Province        NA        NA

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。

# ---- Masking ----
tmax_2018_Jan1_ecu &lt;- 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)
#&gt; 链接到 GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1;sf_use_s2() 为 TRUE
library(terra)
#&gt; terra 1.7.23

# 载入,检查范围
tmax_2018 <- rast("tmax.2018.nc")
ext(tmax_2018)
#&gt; 空间范围 : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)

# 旋转,检查范围
tmax_2018_r <- rotate(tmax_2018)
ext(tmax_2018_r)
#&gt; 空间范围 : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

par(mfrow = c(2, 1))
plot(tmax_2018[[1]])
plot(tmax_2018_r[[1]])

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。<!-- -->


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)

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。<!-- -->

tmax_2018_masked
#&gt; 类型       : SpatRaster 
#&gt; 尺寸  : 13, 34, 365  (nrow, ncol, nlyr)
#&gt; 分辨率  : 0.5, 0.5  (x, y)
#&gt; 空间范围      : -92, -75, -5, 1.5  (xmin, xmax, ymin, ymax)
#&gt; 坐标参考 : 经度/纬度 WGS 84 
#&gt; 来源   : 内存
#&gt; 名称    :    tmax_1,   tmax_2,   tmax_3,   tmax_4,   tmax_5,   tmax_6, ... 
#&gt; 最小值  :  9.768158, 17.85175, 15.38072, 17.51550, 13.29921, 13.41261, ... 
#&gt; 最大值  : 33.109238, 34.98264, 35.73123, 33.20931, 33.43045, 33.31937, ... 
#&gt; 时间        : 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)
#&gt; Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(terra)
#&gt; terra 1.7.23

# load, check extent
tmax_2018 &lt;- rast(&quot;tmax.2018.nc&quot;)
ext(tmax_2018)
#&gt; SpatExtent : 0, 360, -90, 90 (xmin, xmax, ymin, ymax)

# rotate, check extent
tmax_2018_r &lt;- rotate(tmax_2018)
ext(tmax_2018_r)
#&gt; SpatExtent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

par(mfrow = c(2, 1))
plot(tmax_2018[[1]])
plot(tmax_2018_r[[1]])

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。<!-- -->


ecuador_shp &lt;- read_sf(&#39;ecu_shapefiles/ECU_adm1.shp&#39;)
ecuador_shp &lt;- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)

tmax_2018_crop &lt;- crop(tmax_2018_r, ecuador_shp)
plot(tmax_2018_crop[[1]])
plot(ecuador_shp$geometry, add = TRUE)

tmax_2018_masked &lt;- mask(tmax_2018_crop, ecuador_shp)
plot(tmax_2018_masked[[1]])
plot(ecuador_shp$geometry, add = TRUE)

无法遮蔽空间对象,即使将shapefile的CRS转换为与栅格的CRS匹配。<!-- -->

tmax_2018_masked
#&gt; class       : SpatRaster 
#&gt; dimensions  : 13, 34, 365  (nrow, ncol, nlyr)
#&gt; resolution  : 0.5, 0.5  (x, y)
#&gt; extent      : -92, -75, -5, 1.5  (xmin, xmax, ymin, ymax)
#&gt; coord. ref. : lon/lat WGS 84 
#&gt; source(s)   : memory
#&gt; names       :    tmax_1,   tmax_2,   tmax_3,   tmax_4,   tmax_5,   tmax_6, ... 
#&gt; min values  :  9.768158, 17.85175, 15.38072, 17.51550, 13.29921, 13.41261, ... 
#&gt; max values  : 33.109238, 34.98264, 35.73123, 33.20931, 33.43045, 33.31937, ... 
#&gt; 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 &lt;- rast(&quot;tmax.2018.nc&quot;)
tmax_2018 &lt;- rotate(tmax_2018)

ecu &lt;- vect(&#39;ecu_shapefiles/ECU_adm1.shp&#39;)
tmax_2018_ecu &lt;- 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 &lt;- rast(&quot;tmax.2018.nc&quot;)

ecu &lt;- vect(&#39;ecu_shapefiles/ECU_adm1.shp&#39;)
secu &lt;- shift(secu, 360)
tmax_2018_secu &lt;- 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 &lt;- shift(tmax_2018_secu, -360)

huangapple
  • 本文由 发表于 2023年6月15日 02:14:57
  • 转载请务必保留本文链接:https://go.coder-hub.com/76476464.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定