我正在尝试制作一个带有传单的地图,最终将用于显示在地理区域内插的各种估计值.
我为链接数据道歉.我不够聪明,无法弄清楚如何从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")
但得到了错误
Error in rasterFromXYZ(., crs = “+init=epsg:4326 +proj=longlat
+datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0”) : x cell sizes are not regular
这种作品,但没有颜色,所以我看不到交集
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)
Warning message:
In raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) :
‘from’ has no cell values
我还认为如果我使用空间包,可能会计算出交集:
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)
但我明白了
Warning message: In raster::projectRaster(x, raster::projectExtent(x,
crs = sp::CRS(epsg3857))) : ‘from’ has no cell values
那么,实际的问题是:
>这是解决这个问题的合理方法吗?我错过了更好的方法吗?
>我如何从纬度和经度转到带有值的光栅图像以给出阴影,以便我可以看到我正在使用的是什么?
即使我没有直接回答手头的问题,我也会欣喜若狂地获得一些关于在哪里找到类似问题的例子的指导.我还没有找到类似的东西,这可能是一个不知道正确搜索条件的标志.
编辑:
我觉得我可能更接近一步,但仍然不是那么远.
我根据@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对象的文件.
最佳答案 您可以使用表示陆地区域的shapefile屏蔽栅格.您可以使用您的Ohio tract 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)
它并不完全遵循海岸线,但它并不太糟糕.