在水体上移除部分光栅图像

Posted 2023-03-17

【中文标题】在水体上移除部分光栅图像

【英文标题】：Removing part of a raster image over a body of water

【发布时间】：2016-08-24 14:49:16

【问题描述】：

我正在尝试使用leaflet 制作地图，最终将用于显示在地理区域内插值的各种估计值。

对于链接到数据的行为，我深表歉意。我不够聪明，无法弄清楚如何从 Google Docs 导入 .RData 文件。我正在使用的光栅形状可以从https://drive.google.com/file/d/0Bw3leA1Ef_e5RmdKVDNYX0xmS2s/view?usp=sharing 下载 （这是与下面使用的 Coord、rnew 和 cleveland 对象的更新链接）。

我可以得到我想要的地图：

library(tigris) 
library(leaflet)
library(raster)
library(magrittr)
library(dplyr)

load("Coord.Rdata")

rnew <- 
  rasterFromXYZ(
    Coord,
    crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  )

pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(rnew),
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(rnew, colors = pal, opacity = 0.4)

产生

现在，尽管如此美丽，但在伊利湖上显示估算值对我来说没有多大意义。我想删除湖上的光栅图像的部分，或者将其设置为透明......这不会给人留下我正在计算海洋野生动物风险的印象。

我想如果我能找到俄亥俄州的纬度和经度，我也许可以只使用交叉点，从而消除湖上的点。我试着把它放在一起使用

ohio <- tracts(state = "39")

ohio_raster <- 
  attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %>%
  setNames(c("x", "y")) %>%
  mutate(x = sub("^[-]", "", x) %>%
           as.numeric(),
         y = sub("^[+]", "", x) %>%
           as.numeric(),
         layer = 1) %>%
  as.matrix() %>%
  rasterFromXYZ(crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

但出现错误

rasterFromXYZ 中的错误(., crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") : x 像元大小不规则

这种作品，但没有颜色，所以我看不到交叉点

ohio_raster <- 
  raster(ohio, crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

leaflet() %>% addTiles() %>%
  addRasterImage(rnew, colors = pal, opacity = 0.4) %>%
  addRasterImage(ohio_raster, colors = "red", opacity = .4)

警告信息： 在 raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) : 'from' 没有单元格值

（注意：地图已缩小以包括整个俄亥俄州）

我还认为如果我使用空间包可能会计算出交叉点：

spdf <- 
  SpatialPointsDataFrame(Coord[, 1:2], 
                         Coord,
                         proj = CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

ohio_coord <- attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %>%
  setNames(c("x", "y")) %>%
  mutate(x = sub("^[-]", "", x) %>%
           as.numeric(),
         y = sub("^[+]", "", x) %>%
           as.numeric(),
         layer = 1) %>%
  as.matrix() %>%
  as.data.frame()

spdf_ohio <- 
  SpatialPointsDataFrame(ohio_coord[, 1:2], 
                         ohio_coord,
                         proj = CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

spdf_merge <- merge(spdf, spdf_ohio)

leaflet() %>% addTiles() %>%
  addRasterImage(raster(spdf_merge), colors = pal, opacity = 0.4) 

但我明白了

警告消息：在 raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) : 'from' 没有单元格值

所以，实际的问题：

这是解决这个问题的合理方法吗？有没有我想念的更好的方法？ 如何从纬度和经度转换为带有阴影值的光栅图像，以便我可以看到正在处理的内容？

即使我没有得到手头问题的直接答案，我也会欣喜若狂地获得一些关于在哪里可以找到类似问题示例的指导。我还没有找到类似的东西，这可能是不知道正确搜索词的迹象。

编辑：

我觉得我可能更近了一步，但还不是很远。

我根据@IvanSanchez 的推荐下载了克利夫兰地理。只需绘制多边形即可。

这看起来像我想要保留的土地空间。这是一个很好的步骤。

现在我尝试将其与我当前的栅格相交

cleveland <- 
  readOGR(dsn = "cleveland_ohio.osm2pgsql-shapefiles",
          layer = "cleveland_ohio_osm_polygon",
          p4s = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

leaflet() %>% addTiles() %>%
  addRasterImage(raster::intersect(rnew, cleveland), 
                 colors = pal, 
                 opacity = 0.4)

但我似乎拥有联合，而不是区域的交叉点。或者不同的东西。

我已将上面的链接更新为包含Coord、rnew 和cleveland 对象的文件。

【参考方案1】：

您可以使用代表陆地区域的 shapefile 来屏蔽栅格。您可以使用俄亥俄州的 shapefile，但 Natural Earth 的 lakes data 可能会简化事情。例如：

library(leaflet)
library(raster)
library(magrittr)
library(rgdal)
library(rgeos)

load('Coord.Rdata')  
rnew <- rasterFromXYZ(Coord, crs='+init=epsg:4326')
pal <- colorNumeric(c('#0C2C84', '#41B6C4', '#FFFFCC'), values(rnew),
                    na.color = 'transparent')

# Download and extract the Natural Earth data
download.file(
  file.path('http://www.naturalearthdata.com/http/',
            'www.naturalearthdata.com/download/10m/physical/ne_10m_lakes.zip',
            f <- tempfile())
unzip(f, exdir=tempdir())

# Read in the lakes shapefile    
lakes <- readOGR(tempdir(), 'ne_10m_lakes')

# Create a "land" shapefile by taking the difference between the 
# raster's extent and the lakes polys
land <- gDifference(as(extent(rnew), 'SpatialPolygons'), lakes)

# Mask rnew
rnew2 <- mask(rnew, land)

leaflet() %>% addTiles() %>%
  addRasterImage(rnew2, colors = pal, opacity = 0.4)

它并不完全遵循海岸线，但也不算太糟糕。

【参考方案2】：

我建议您从MapZen metro extract for Cleveland 下载该地区的地理数据。这将为您提供一些用于陆地和水域的漂亮多边形 - 然后您可以将您的形状与陆地多边形相交（或从您的形状中减去水多边形）。

【讨论】：

在尝试了您的建议后，我做了一些修改。我得到了克利夫兰的形状，但仍然不能只拉出十字路口。看起来raster::intersection 是正确的做法，但显然我还是什么都没看到。