第 6 章 图例

6.1 6行代码解决ArcGIS栅格数据配色

如果你是ArcGIS的用户,你肯定被ArcGIS的栅格图像的配色Classified Symbology折磨过。

这是一个典型的ArcGIS style的图例,虽然奇丑,但是依然需要耗费你巨大的精力做出这个丑东西。

工作量在于:

  1. breaks不能自动设定,必须手动一个个敲进去。分级发生改变之后,之前敲的breaks又全部乱掉。
  2. labels同breaks一样,每次微调,必须老老实实breaks和labels全部手动敲一遍。
  3. 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
val_max = 9999
brks_land = c(0, 2, 5, 10, 50, 100, 1000, val_max) # 陆地高程breaks
brks_sea  = c(-val_max, seq(-60, -10, 10), -5, 0)  # 海洋高程breaks
brks = c(brks_sea, brks_land) %>% unique()

# brks2用来制作legend,get_colorkey2碰到Inf或-Inf,才会显示三角尖尖(ncl style)。
brks2 = brks
brks2[brks2 ==  val_max] = Inf
brks2[brks2 == -val_max] = -Inf 

6.1.2 colormap

选色是个困难的技术活。有利器为你助力:(1) GrADS-NCL_colormap_V1.3,该软件可以自动从别人图件的legend中copy图例颜色;

  1. R包rcolors为你提供270种colormaps。
read_color <- function(file){
    tryCatch({
        d = read.table(file, skip = 1) %>% set_names(c("r", 
            "g", "b"))
        max_value = max(sapply(d, max))
        max_value = ifelse(max_value <= 1, 1, 255)
        colors = with(d, rgb(r, g, b, maxColorValue = max_value))
        colors
    }, error = function(e) {
        message(sprintf("[%s]: %s", basename(file), e$message))
    })
}

resample_color <- function(cols, n = NULL, show = FALSE) 
{
    # cols = rcolors[[name]]
    if (is.null(n)) 
        n = length(cols)
    cols = colorRampPalette(cols)(n)
    if (show) 
        show_col(cols)
    cols
}
## 2. cols
blues = RColorBrewer::brewer.pal(9, "Blues")[-(1:2)] %>% rev()# %>% show_cols()
cols_sea  = resample_color(blues, length(brks_sea)-1) 

file_landcolor = "8colors.rgb" # 此为利用`GrADS-NCL_colormap_V1.3`copy来的colors
cols_land = rcolors::read_color(file_landcolor)[-1]
#> Warning in file(file, "rt"): cannot open file '8colors.rgb': No such file or
#> directory
#> [8colors.rgb]: cannot open the connection
cols = c(cols_sea, cols_land)

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
write_sp2rgb <- function(grid, brks, cols, file = "sp_rgb.tif"){
    colormap = col2rgb(cols) %>% t() %>% as.data.table()
    ind   = findInterval(grid$band1, brks) # %>% summary()
    d_rgb = colormap[ind, ]
    
    grid_rgb   = grid
    grid_rgb@data <- d_rgb
    writeGDAL(grid_rgb, file, options = c("COMPRESS=DEFLATE")) # 
}

TIFF导入导出,两行搞定,自由任性。

# file_dem = "G:/ArcGIS_common/Projects/raster/dem_srtm_gz2"
grid <- readGDAL("data-raw/dem_pearl.tif" )
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
key = latticeMap::get_colorkey(brks2, cols, is_factor = TRUE)
g <- draw.colorkey(key)
grid::grid.draw(g)

write_fig(g, "images/dem_legend.jpg", 9.2, 0.75, show = FALSE)

6.1.5 ArcGIS中拼接

最终拼图如下所示,中间过程不再赘述。