倍频程:盒形和晶须图，未从GeoTIFF传播 [英] Octave: Box and whisker plot without spread from GeoTIFF

问题描述

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

使用八度音阶，我想创建30个水平的箱形图和晶须图，并且没有来自30个不同的GeoTIFF的扩散(x轴).这是我希望情节的样子的草图:

Using Octave I want to create 30 horizontal box and whisker plots without spread (x-axis) from 30 different GeoTIFF's. This is a sketch of how I would like the plot to look like:

理想地，对我来说最好的解决方案是八度代码(工作流程)，该代码使我可以在一个目录中放置多个GeoTIFF，然后单击一下即可一次为所有GeotIFF创建一个框和晶须图-就像上面的草图一样

Ideally the best solution for me would be an Octave code (workflow) that would allow me to place multiple GeoTIFFs in one directory and then with one click create a box and whisker plot for all GeotIFFs at once - just like the sketch above.

可以在此处下载具有3个GeoTIFF的GeoTIFF示例.一个>.该文件在QGIS中如下所示:

A GeoTIFF-sample with 3 GeoTIFF's can be downloaded here. The file looks like this in QGIS:

它保留了波段1上的高程值(每个箱形图和晶须图应基于该高程值，并且没有数据值(-999)，该图上应排除无数据值.

It holds elevation values on band 1 (the ones that each box and whisker plot should be based on, and no data values (-999), the no-data values should be excluded from the plot.

现在这就是我得到的:

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

2. 要加载它并创建箱形图，我目前正在使用此代码(感谢@Andy和@Cris Luengo):

1. Using img = imread ("filname.tif") gets the file into Octave. Using hist (img(:), 200); shows that all cells are concentrated around 65300. imagesc (img, [65100 65600]) follwed by colorbar displays the image extent but's it's clear that this way simply doesn't import the real cell values. I can't find a working solution to import GeoTIFF's with cell values, therefor my current work-around is exporting the GeoTIFF from QGIS with gdal_translate -of aaigrid which creates a .asc-file that I manually edit to remove header rows, rename to .csv and load into Octave. That .csv can be found here.

2. To load it and create a box plot I'm currently using this code (thanks to @Andy and @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轴(从200-450开始，而不是相反).这是由view([-90 90])引起的，我使用该view([-90 90])来使箱形图水平放置，而不是出于布局原因而需要垂直放置.

The results is pretty close but I'm still failing to: A) load GeoTIFF's directly from a folder. If this is not possible I'm gonna have to modify the code to load all *.csv in a directory to the same box plot and label each plot by filename (which I'm unsure how to accomplish. B) to get the x-axis reversed (going from 200-450, not the other way around). This is caused by the view([-90 90]) that I use to make the box plot horizontal instead of vertical which is needed for layout reasons.

有人对如何解决最近的调整有任何想法吗?

Anyone with any ideas on how to resolve the last adjustments?

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

我有30个GeoTIFF，其中包含来自视域分析的结果，每2x2米见方，有一个值告诉我在从视域点可见之前建筑物的高度(以米为单位).结果覆盖了整个斯德哥尔摩市，但是上述30个GeoTIFF是计划进行新开发的区域的较小片段.结果有助于规划人员了解新开发如何影响这30个地方(这对文化遗产管理非常重要).

I have 30 GeoTIFFs containing results from a viewshed analysis, for every 2x2 meter square there is a value the tells me how high a building can be (in meters) before it's visible from the viewshed point. The results cover the whole city of Stockholm but the above mentioned 30 GeoTIFFs are smaller clips of an area where new development is planned. The results help planners to understand how new development might effect each of the 30 places (that are important for cultural heritage management).

作为更大的PDF报告的一部分(这些结果可以用不同比例的不同地图可视化)，我试图制作箱形图和晶须图(作为对地图的补充)，从而为读者提供了一个概述根据30个视域(GeoTIFF)的每个结果(30个位置的每个盒子和一个晶须)，规划的开发区域还剩下多少空间.以下是报告中的地图外观的示例:

As part of a bigger PDF-report (where these results are visualized with different maps in different scales) I'm trying to produce a box and whisker plot (as a compliment to the maps) the gives the reader an overview over how much space is there is left at the planned development area, based on each of the 30 viewshed (GeoTIFF) results (one box and whisker for each of the 30 locations). Below is an example of how a map in the report can look like:

推荐答案

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

Does not directly read GeoTIFF but calls gdal_translate under the hood. Just place all your .tif in the same directory. Make sure gdal_translate is in your 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")

给予

这篇关于倍频程:盒形和晶须图，未从GeoTIFF传播的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！