八度：没有从 GeoTIFF 传播的盒须图

Question

---- 可以在下面的第 3 点中找到到目前为止我所得到的更新以及需要解决的问题----

我想使用 Octave 从 30 个不同的 GeoTIFF 中创建 30 个没有展开的水平盒须图 (x-axis)。这是我希望情节看起来像的草图：

理想情况下，对我来说最好的解决方案是 Octave 代码（工作流程），它允许我将多个 GeoTIFF 放在一个目录中，然后单击一次为所有 GeotIFF 创建一个盒须图 - 就像上图。

可以下载具有 3 个 GeoTIFF 的 GeoTIFF-sample here。该文件在 QGIS 中看起来像这样：

它保存波段 1 的高程值（每个盒须图应基于这些值，并且没有数据值 (-999)，no-data 值应从图中排除。

现在这是我得到的：

使用 img = imread ("filname.tif") 将文件放入 Octave。使用 hist (img(:), 200); 显示所有单元格都集中在 65300 左右。imagesc (img, [65100 65600]) 后跟 colorbar 显示图像范围，但很明显这种方式根本不会导入真实的单元格值。我找不到使用单元格值导入 GeoTIFF 的工作解决方案，因此我当前的 work-around 正在使用 gdal_translate -of aaigrid 从 QGIS 导出 GeoTIFF，这会创建一个 .asc-file 我手动编辑以删除header 行，重命名为 .csv 并加载到 Octave。那个.csv可以找到here.

要加载它并创建一个箱线图，我目前正在使用此代码（感谢@Andy 和@Cris Luengo）：

pkg load statistics

s = urlread ("https://drive.google.com/uc?export=download&id=1RzJ-EO0OXgfMmMRG8wiCBz-51RcwSM5h");
o = str2double (strsplit (s, ";"));
o(isnan (o)) = [];

boxplot (o)
set(gca,"xtick",[])
view([-90 90])

print out.png

结果非常接近，但我仍然无法：A) 直接从文件夹加载 GeoTIFF。如果这不可能，我将不得不修改代码以将目录中的所有 *.csv 加载到同一个箱线图，并按文件名标记每个图（我不确定如何完成。B) 使 x-axis 反转（从 200-450，而不是相反）。这是由我用来使箱线图水平而不是垂直的 view([-90 90]) 引起的，这是布局原因所需要的。

有人对如何解决最后的调整有任何想法吗？

----背景信息----

我有 30 个包含视域分析结果的 GeoTIFF，每 2x2 平方米都有一个值告诉我建筑物在从视域点可见之前可以有多高（以米为单位）。结果涵盖了整个斯德哥尔摩市，但上述 30 个 GeoTIFF 是规划新开发区域的较小片段。结果有助于规划人员了解新开发项目可能如何影响 30 个地方（对文化遗产管理很重要）。

作为更大 PDF-report 的一部分（这些结果用不同比例的不同地图可视化）我正在尝试制作一个盒须图（作为对地图的补充）给出了reader 根据 30 个视域 (GeoTIFF) 结果中的每一个（30 个位置中的每一个位置都有一个方框和胡须），space 计划开发区域还剩多少的概览。以下是报告中地图的外观示例：

Answer 1

不直接读取 GeoTIFF，而是在后台调用 gdal_translate。只需将所有 .tif 文件放在同一个目录中。确保 gdal_translate 在你的路径中：

pkg load statistics

clear all;
fns = glob ("*.tif");
for k=1:numel (fns)

  ofn = tmpnam;
  cmd = sprintf ('gdal_translate -of aaigrid "%s" "%s"', fns{k}, ofn);
  [s, out] = system (cmd);
  if (s != 0)
    error ('calling gdal_translate failed with "%s"', out);
  endif
  fid = fopen (ofn, "r");
  # read 6 headerlines
  hdr = [];
  for i=1:6
    s = strsplit (fgetl (fid), " ");
    hdr.(s{1}) = str2double (s{2});
  endfor
  d = dlmread (fid);

  # check size against header
  assert (size (d), [hdr.nrows hdr.ncols])

  # set nodata to NA
  d (d == hdr.NODATA_value) = NA;

  raw{k} = d;

  # create copy with existing values
  raw_v{k} = d(! isna (d));

  fclose (fid);

endfor

## generate plot
boxplot (raw_v)
set (gca, "xtick", 1:numel(fns),
          "xticklabel", strrep (fns, ".tif", ""));
view ([-90 90])
zoom (0.95)

print ("out.png")

给予

八度：没有从 GeoTIFF 传播的盒须图

Octave: Box and whisker plot without spread from GeoTIFF

octave

geotiff