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

Question

---- 更新我目前得到的内容，还有待解决的问题可以在下面的第 3 点找到----

使用 Octave，我想从 30 个不同的 GeoTIFF 中创建 30 个没有散布（x 轴）的水平箱线图和晶须图。这是我希望情节看起来如何的草图：

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

可以下载带有 3 个 GeoTIFF 的 GeoTIFF 样本here。该文件在 QGIS 中如下所示：

它保存带 1 上的高程值（每个箱线图应该基于的高度值，并且没有数据值 (-999)，应该从图中排除无数据值。

现在这是我得到的：

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

2. 要加载它并创建一个箱形图，我目前正在使用此代码（感谢 @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
   

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

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

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

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

作为更大的 PDF 报告的一部分（这些结果通过不同比例的不同地图进行可视化），我试图制作一个箱线图（作为对地图的补充），让读者对根据 30 个视域 (GeoTIFF) 结果中的每一个（30 个位置中的每一个都有一个方框和胡须），规划的开发区域还剩下多少空间。以下是报告中地图的外观示例：

Answer 1

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

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")

给予