八度:没有从GeoTIFF传播的盒子和胡须情节

Question

使用八度，我想从30个不同的GeoTIFF中创建30个水平框和胡须情节，而不需要展开(x轴)。这是我希望这个情节看起来如何的草图：

​

理想情况下，对我来说最好的解决方案是一个Octave代码(工作流)，它允许我将多个GeoTIFFs放在一个目录中，然后单击一次，为所有GeotIFF创建一个框和晶须图--就像上面的草图一样。

一个GeoTIFF-有3个GeoTIFF的样本可以下载这里。该文件在QGIS中如下所示：

​

它保存在波段1上的高程值(每个框和晶须图应该基于的那些值，而没有数据值(-999)，不应该从图中排除非数据值。

现在，这就是我得到的：

1. 使用img = imread ("filname.tif")将文件输入到Octave中。使用hist (img(:), 200);显示，所有细胞集中在65300左右。imagesc (img, [65100 65600])由colorbar折叠，显示图像范围，但很明显，这种方式根本不导入真正的单元格值。我找不到一个工作解决方案来导入带有单元格值的GeoTIFF，因此我当前的工作是用gdal_translate -of aaigrid从QGIS导出gdal_translate -of aaigrid，它创建了一个..asc文件，我手动编辑该文件以删除标题行，重命名为.csv并加载到Octave中。.csv可以被找到这里。
2. 要加载它并创建一个方框图，我目前正在使用以下代码(感谢@Andy和@Cris Luengo)： pkg负载统计信息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"，[])视图(- 90)打印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在您的路径上：

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

给出

​

​