第 6 章 图例
6.1 6行代码解决ArcGIS栅格数据配色
如果你是ArcGIS的用户,你肯定被ArcGIS的栅格图像的配色Classified Symbology折磨过。
这是一个典型的ArcGIS style的图例,虽然奇丑,但是依然需要耗费你巨大的精力做出这个丑东西。
工作量在于:
- breaks不能自动设定,必须手动一个个敲进去。分级发生改变之后,之前敲的breaks又全部乱掉。
- labels同breaks一样,每次微调,必须老老实实breaks和labels全部手动敲一遍。
- ArcGIS的colors更是死板,template有限且固定。自动生成的颜色一般很难满足要求。这样的话,每个level你都要手动选择颜色,每次只能瞪大眼睛试试这个颜色怎么样。
ArcGIS虽然有python arcpy
,但google搜遍,这么简单的需求它却无法实现。但是ArcGIS其他方面的灵活,又是MATLAB、R无法替代的,因此又爱又恨又离不开它。
本文介绍一种非常简易的方法,配合R语言进行数据预处理,上述过程中简单到不可思议,6行代码解决栅格图像的配色、3行代码造ncl风格图例。
众所周知,TIFF可以包含大地坐标系统信息,TIFF同时又可以显示RGB 3波段彩色图像。据此,我们可以把配好颜色的栅格保存成具有空间位置的RBG类型的TIFF。先根据栅格数据的values,在不同的breaks区间中填充不同的颜色,然后在把RGB值写入到原栅格位置。最后再徒手造一个legend,塞给ArcGIS。这样就实现了根据breaks填充colors,又实现了保留空间位置信息,因此导入ArcGIS无需校准空间位置。
6.1.1 造数据
以珠江三角洲高程专题图为例
巧妇难为无米之炊,先造breaks和colors:
## brks in raster::plot can not have Inf values
= 9999
val_max = c(0, 2, 5, 10, 50, 100, 1000, val_max) # 陆地高程breaks
brks_land = c(-val_max, seq(-60, -10, 10), -5, 0) # 海洋高程breaks
brks_sea = c(brks_sea, brks_land) %>% unique()
brks
# brks2用来制作legend,get_colorkey2碰到Inf或-Inf,才会显示三角尖尖(ncl style)。
= brks
brks2 == val_max] = Inf
brks2[brks2 == -val_max] = -Inf brks2[brks2
6.1.2 colormap
选色是个困难的技术活。有利器为你助力:(1) GrADS-NCL_colormap_V1.3,该软件可以自动从别人图件的legend中copy图例颜色;
- R包rcolors为你提供270种colormaps。
<- function(file){
read_color tryCatch({
= read.table(file, skip = 1) %>% set_names(c("r",
d "g", "b"))
= max(sapply(d, max))
max_value = ifelse(max_value <= 1, 1, 255)
max_value = with(d, rgb(r, g, b, maxColorValue = max_value))
colors
colorserror = function(e) {
}, message(sprintf("[%s]: %s", basename(file), e$message))
})
}
<- function(cols, n = NULL, show = FALSE)
resample_color
{# cols = rcolors[[name]]
if (is.null(n))
= length(cols)
n = colorRampPalette(cols)(n)
cols if (show)
show_col(cols)
cols }
## 2. cols
= RColorBrewer::brewer.pal(9, "Blues")[-(1:2)] %>% rev()# %>% show_cols()
blues = resample_color(blues, length(brks_sea)-1)
cols_sea
= "8colors.rgb" # 此为利用`GrADS-NCL_colormap_V1.3`copy来的colors
file_landcolor = rcolors::read_color(file_landcolor)[-1]
cols_land #> Warning in file(file, "rt"): cannot open file '8colors.rgb': No such file or
#> directory
#> [8colors.rgb]: cannot open the connection
= c(cols_sea, cols_land) cols
6.1.3 栅格配色
#' write_sp2rgb
#'
#' @param grid SpatialGridPixels or SpatialGridDataFrame object
#' @param brks breaks
#' @param cols colors
#' @param file output filename of tif file
#'
#' @examples
#' write_sp2rgb(grid, brks, cols, file = "dem_pearl_rgb.tif")
#' @export
<- function(grid, brks, cols, file = "sp_rgb.tif"){
write_sp2rgb = col2rgb(cols) %>% t() %>% as.data.table()
colormap = findInterval(grid$band1, brks) # %>% summary()
ind = colormap[ind, ]
d_rgb
= grid
grid_rgb @data <- d_rgb
grid_rgbwriteGDAL(grid_rgb, file, options = c("COMPRESS=DEFLATE")) #
}
TIFF导入导出,两行搞定,自由任性。
# file_dem = "G:/ArcGIS_common/Projects/raster/dem_srtm_gz2"
<- readGDAL("data-raw/dem_pearl.tif" )
grid write_sp2rgb(grid, brks, cols, file = "dem_pearl_rgb.tif")
6.1.4 图例
关于labels,它的主要用途是制作图例。我们完全可以采用其他软件制作图例,保存成jpg。然后再导回ArcGIS即可。这样就避开了在ArcGIS中设置。Ipaper为你提供了对应的函数get_colorkey2
,它可以制作ncl style的colorbar。
library(Ipaper)
# 若无get_colorkey2函数,记得先更新下Ipaper
= latticeMap::get_colorkey(brks2, cols, is_factor = TRUE)
key <- draw.colorkey(key)
g ::grid.draw(g) grid
write_fig(g, "images/dem_legend.jpg", 9.2, 0.75, show = FALSE)
6.1.5 ArcGIS中拼接
最终拼图如下所示,中间过程不再赘述。