在Octave中导入和编写GeoTIFF

Question

我在办公室使用MATLAB，在家里用Octave。虽然它们非常相似，但我试图做一些看起来很简单和显而易见的事情，但却发现真的很烦人。我找不到如何在Octave中导入TIFF图像的方法。我知道MATLAB的geotiffread函数不存在，但我认为会有另一种方法。

在某些情况下，我可以跳过导入它们，因为我可以使用imread函数进行操作，但第二个问题是我找不到一种编写带地理参考的TIFF文件的方法（在MATLAB中，我通常在geotiffwrite中调用geotiffinfo输入）。我的TIFF文件通常是8位无符号整数或32位有符号整数。我希望有人能提出解决这个问题的方法。我还看到了此主题，但不确定是否可以在Octave中使用Ashish提出的代码。

- umbe1987

---

看起来未维护的映射包有你正在寻找的geotiffread函数。 - juliohm

我认为映射包还没有完善，它仍然缺少一个函数。 - umbe1987

请参见下面为其编写C++包装器的答案。 - makarand kulkarni

Answer 1

你可以查看Octave中的映射库。你也可以使用栅格函数来处理GeoTiffs。

示例：

pkg load mapping
filename=”C:\\sl\SDK\\DTED\\n45_w122_1arc_v2.tif”
rasterinfo (filename)
rasterdraw (filename)

- Jeremy Sarao

Answer 2

简化版：Octave默认情况下无法处理有符号32位图像，但这并非不可能。由于Octave是免费软件，它的功能取决于用户愿意花费时间或金钱来实现哪些功能。从版本3.8.1开始，Octave使用GraphicsMagick或ImageMagick来处理图像的读写，但存在一些问题，如精度受限和只能写无符号整数等。对于geotiff格式的处理也存在一些问题。

详细版：短答案是你不能在Octave中直接处理有符号32位图像。但这并不是因为不可能做到，而是因为还没有人费心去实现它。作为一款免费软件，Octave的功能取决于其用户愿意花费时间或金钱来实现哪些功能。从版本3.8.1开始，Octave使用GraphicsMagick或ImageMagick来处理图像的读写。但这引入了一些问题。其中一个问题是，你的精度受制于你构建GraphicsMagick时的量子深度选项。此外，你只能写入无符号整数。希望这种情况在未来会有所改变，但由于没有多少用户需要它，所以一直保持这样的状态。对于geotiff格式的处理也存在一些问题。

只要你了解C++，就可以自己编写这些函数。这应该不会太难，因为已经有了一个名为libgeotiff的C库。你只需要编写一个Octave oct函数的包装器（当然，如果你不懂C或C++，那么这个“只”就变成了很多工作）。

- carandraug

Answer 3

这是需要编译的示例oct文件代码。我参考了https://gerasimosmichalitsianos.wordpress.com/2018/01/08/178/。

#include <octave/oct.h>
#include "iostream"
#include "fstream"
#include "string"
#include "cstdlib"
#include <cstdio>

#include "gdal_priv.h"
#include "cpl_conv.h"
#include "limits.h"
#include "stdlib.h"
using namespace std;
typedef std::string String;

DEFUN_DLD (test1, args, , "write geotiff")
{
  NDArray maindata = args(0).array_value ();
  const dim_vector dims = maindata.dims ();
  int i,j,nrows,ncols;
    nrows=dims(0);
    ncols=dims(1);  
        //octave_stdout << maindata(i,0);
    NDArray transform1 = args(1).array_value ();
  double*  transform = (double*) CPLMalloc(sizeof(double)*6);
  float*  rowBuff = (float*) CPLMalloc(sizeof(float)*ncols);
  //GDT_Float32 *rowBuff = CPLMalloc(sizeof(GDT_Float32)*ncols);
  String tiffname;
  tiffname  = "nameoftiff2.tif";
  
  cout<<"The transformation matrix is";
  for (i=0; i<6; i++)
  {
      transform[i]=transform1(i);
      cout<<transform[i]<<" ";
  }
     
  GDALAllRegister();
  CPLPushErrorHandler(CPLQuietErrorHandler);
  GDALDataset *geotiffDataset;
  GDALDriver *driverGeotiff;
  GDALRasterBand *geotiffBand;
  
  OGRSpatialReference oSRS;
  char **papszOptions = NULL;
  char *pszWKT = NULL;
 
    oSRS.SetWellKnownGeogCS( "WGS84" );
    oSRS.exportToWkt( &pszWKT );
    
  driverGeotiff = GetGDALDriverManager()->GetDriverByName("GTiff");  
  geotiffDataset = (GDALDataset *) driverGeotiff->Create(tiffname.c_str(),ncols,nrows,1,GDT_Float32,NULL);
  geotiffDataset->SetGeoTransform(transform);
  geotiffDataset->SetProjection(pszWKT);
  //CPLFree( pszSRS_WKT );
  cout<<" \n Number of rows and columns in array are: \n";
  cout<<nrows<<" "<<ncols<<"\n";
  for (i=0; i<nrows; i++)
    {
        for (j=0; j <ncols; j++)
            rowBuff[j]=maindata(i,j);
        //cout<<rowBuff[0]<<"\n";
        geotiffDataset->GetRasterBand(1)->RasterIO(GF_Write,0,i,ncols,1,rowBuff,ncols,1,GDT_Float32,0,0);
    }
  
  
  GDALClose(geotiffDataset) ;
  CPLFree(transform);
  CPLFree(rowBuff);
  CPLFree(pszWKT);
  GDALDestroyDriverManager();
  
  return octave_value_list();
}

可以使用以下方式进行编译和运行

mkoctfile -lgdal test1.cc
aa=rand(50,53);
b=[60,1,0,40,0,-1];
test1(aa,b);

- makarand kulkarni