python-GDAL空间分析编程的整理

记得前段时间，总结了基于python的GIS空间分析类库，大大小小有20多个，司职于不同的领域，当然也有一定的重复，不过对于空间分析部分，总结下来最为主要的还是python-GDAL以及矩阵分析类库numpy，ArcGIS 最新版本已经将numpy囊括到它的geoprocessing里面去了，看来用numpy做空间分析是大势所趋。因为一直有用基于python写Multi-agents based land use modelling 的需求， 所以在空间分析方面有了一点小小的体会和经验，不成系统，希望大家多多探讨交流。

当然，如果要做空间分析，首先要将数据读入，gdal是必须的，要分析才牵扯到numpy做栅格分析，shapely做矢量分析。记得以前还列出了scipy等，按照长尾理论，这些类库用的机会是很少的，所以我在想，如果要开发一整套基于python的空间类库，是不是有必要引入这个庞大的统计分析模块。自己读过UCSB一些牛人写的程序，他们甚至连数据读入和输出都是自己来写，我方是觉得这似乎没有必要，基于GDAL在GIS业界的声誉和强大功能，引入GDAL是不会错的。

要做空间分析，就牵扯到数据输入、处理，输出三个部分。关于数据输入，我情钟与GeoTiff，虽然大牛们用asc格式（AAGrid）比较多，但是GeoTiff的优势是显而易见。至少在gdal的GTiff driver的create和createCopy方法都会支持。

1
2
3
4
5
6
7
8
9

    format = "HFA"
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) \
       and metadata[gdal.DCAP_CREATE] == 'YES':
        print 'Driver %s supports Create() method.' % format
    if metadata.has_key(gdal.DCAP_CREATECOPY) \
       and metadata[gdal.DCAP_CREATECOPY] == 'YES':
        print 'Driver %s supports CreateCopy() method.' % format

其次，在数据读入这里，我建议自己先确定整个分析研究区的情况，比如大小，投影，转换，控制点什么的，这些metadata数据最好是import一个基本数据，比如import base year的土地利用情况。其实对于这点，大家可以思量以前用ArcMap加载数据的情况，如果要设置一个layer的相关信息，包括投影信息，extent信息等，我们都会选择导入项目中和该数据具有相同metadata的数据。通过以上分析，可以看出，如果要建立一个基本研究区，可以通过以下两种方法：

• 自己制作数据的metadata，相对而言，比较麻烦，
• 导入相关信息的metadata，方便，不过需要有相关信息，如果没有，只有用前者。

例如,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

    def init_from_name(self, filename):
        """
        Initialize file_info from filename
        filename -- Name of file to read.
        Returns 1 on success or 0 if the file can't be opened.
        """
        fh = gdal.Open( filename )
        if fh is None:
            return 0
        self.filename = filename
        self.bands = fh.RasterCount
        self.xsize = fh.RasterXSize
        self.ysize = fh.RasterYSize
        self.band_type = fh.GetRasterBand(1).DataType
        self.projection = fh.GetProjection()
        self.geotransform = fh.GetGeoTransform()
        self.ulx = self.geotransform[0]
        self.uly = self.geotransform[3]
        self.lrx = self.ulx + self.geotransform[1] * self.xsize
        self.lry = self.uly + self.geotransform[5] * self.ysize
 
        ct = fh.GetRasterBand(1).GetRasterColorTable()
        if ct is not None:
            self.ct = ct.Clone()
        else:
            self.ct = None
        return 1

这样就有了一个基本研究区概况的类，这样我们做分析的时候，就不用不停的去获取metadata信息，另一个需要注意的是关于nodatavalue，这个以前在用erdas和arcgis做分析的时候并没有太多注意，不过如果要用python编程，nodatavalue的值是个很重要的部分，搞不好，你在做slope分析或者做hillshade的时候，把不需要分析的值加以分析，产生错误的结果。一般情况下，如果是和邻域无关的空间分析还好点，如果和邻域有关，那结果肯定是有问题的。比如说，你要做一个moore准则的CA，如果你不设定nodatavlue那么，这样一来所有做为nodata的部分全部都死亡掉了，然后进入了整个loop。

对于nodatavalue，我个人的看法是，首先要有一个整个研究区的nodatavlue，这样方便我们输出结果，其次很多情况下，栅格数据不光包括DEM数据，（float数据，一般在-8000～8000之前），slope数据（90以下），AHP法和delphi法得到的打分数据（基于打分标准，可能有float 0~1, float 0~100, byte 0~255， int 0～1000 or 10000）,光栅值数据（如果是一个band，那么就是byte 0～255）. 这样看来，还是找不到一个合理的nodatavalue，因为大家的nadata的range实在是太多了，比如说用 数据类型的最大或者最小值做nodatavalue的。
GDAL data type minimum maximum
Byte 0 255
UInt16 0 65,535
Int16, CInt16 -32,768 32,767
UInt32 0 4,294,967,295
Int32, CInt32 -2,147,483,648 2,147,483,647
Float32, CFloat32 -3.4E38 3.4E38
Float64, CFloat64 -1.79E308 1.79E308
比如，如果是Byte，就用 0 或者 255， Int16 就用 -32768 或者32767.针对此，我看了gdal的python sample和perrygeo自己推荐的一些实例，（一句题外话，这两个项目的python代码推荐有空大家读一读，还蛮有启发的）,最好的方法是Band.GetNoDataValue()，然后输出的时候再Set.NoDataValue()。

刚才扯到了栅格数据的range问题，当然具有具体地理意义的栅格数据的range我没法去设定，但是关于打分数据，我觉得我们可以规定一个范围，比如说0～1，然后数据类型用float，切记，这只是个打分数据，数据输出的时候，可能还要进行分数在加以reclassify。其实用栅格进行处理打分只是一个最基本的问题，相关问题还有进行标准化，离散化，空间运算等。这些时候，就很难控制一个值域。所以我觉得，不能去控制输入数据的值域，我们能做的，就是在读入数据后进行必要的标准化运算。

最后说说输出，PIL，Pylab，gdal，mapnik都有产生图片的方法。总而言之，如果你还要进行地理分析，或者不想引入其他类库，gdal，PIL就够用了。

1
2
3
4
5
6
7
8
9
10
11
12

    # Output the Array
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    dst_ds = driver.Create( outfile, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32 )
    if ds.GetGeoTransform():
        dst_ds.SetGeoTransform(ds.GetGeoTransform())
    if ds.GetMetadata():
        dst_ds.SetMetadata(ds.GetMetadata())
    if ds.GetProjection():
        dst_ds.SetProjection(ds.GetProjection()) 
    dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )   
    dst_ds.GetRasterBand(1).WriteArray( outArray )

Pylab是用来进行数值分析，和统计分析时使用的。

1
2
3

pylab.imshow(result)
pylab.savefig("/XXX.png")
pylab.draw()

如果你要进行美化渲染并以一种很好的形式表现，mapnik是个不错的选择，现在mapnik已经更新到0.6.0了，变化还是很大的，有空给大家介绍介绍。

最后说说架构，建议大家可以读读几个项目的源代码，比如cage，urbansim，gdal python sample。如果你只是先实现一些小功能，那么看看perrygeo的python tools还是挺有帮助的。 就先说这么多，有空再和大家聊。

Related posts:

• 人人都说我丑，其实我只是美得不明显。
• One-day-trip in python for GIS
• 基于python的栅格矢量转换
• 关于MapServer的Output Format
• GDAL/OGR 1.7.0 RELEASED