SuperMapPy 批量拼接 GeoTiff影像

影像拼接工具使用说明

 

一、影像像素位深检查

 

1.采用开源库GDAL的gdalinfo.exe读取GeoTiff文件的信息，如类型、投影，范围等。

2.采用DOS脚本遍历所有GeoTiff文件，输出各个GeoTiff文件的信息到日志。

3.统计GeoTiff文件个数，Type=Byte的文件的个数。数目一致表明全部GeoTiff像素类型一致。

 

运行示例：

for /f %i in ('dir /b /s d:\tifdata\*.tif ') do @(

    e:\gdal1.4.2\bin\gdalinfo.exe %i >> e:\tifstats.log

    )

 

 

二、影像拼接处理

采用SuperMapPy脚本处理影像拼接。先安装Python2.7.3和Deskpro环境。脚本如下：

# coding: GB2312

 

#===================================================

#影像成批导入SuperMap UDB格式工具.

#基本流程:

#1、扫描目录，或者根据经纬度按照块的高宽生成文件列表。

#2、扫描存在的文件，获取最大的坐标范围。

#3、扫描存在的文件，获取像素格式。

#4、遍历文件列表，对于存在的文件追加到打开的数据库UDB中。

#5、创建金字塔索引，以加快显示速度。（可选的过程）

#===================================================

 

import sys

import string

import re

import os

import time

 

def getType(ext):

    if ext.lower() == 'tif':

        return 'fileTIF'

    elif ext.lower() == 'img':

        return 'fileIMG'

 

 

 

#匹配正则表达式，符合条件的append到datafiles，用于追加   

def walkPath(type, path):

    datafiles = []

    reMatch = '[\d\D]*.tif$'

    if type=='img':

        reMath = '[\d\D]*.img$'

   

    for root, dirs, files in os.walk(path):

        for file in files:

            if (re.match(reMatch,file)):

                datafiles.append(os.path.join(root, file))

    print len(datafiles)               

    return datafiles

 

def calcDatasetInfo(type, datafiles):

    L=[]       

    left=[]

    top=[]

    right=[]

    bottom=[]

    ratiox=[]

    ratioy=[]

 

    #获取每个影像文件的左右地理范围，保存到数组

    for file in datafiles:

        L= smu.GetImageGeoRef(type,file)

        print L

        l=float(L[0][0])

        t=float(L[0][1])

        r=float(L[0][2])

        b=float(L[0][3])

        w=int(L[1][0])

        h=int(L[1][1])

        x=(r-l)/w

        y=(t-b)/h

       

        left.append(l)

        right.append(r)

        top.append(t)

        bottom.append(b)

        ratiox.append(x)

        ratioy.append(y)

 

        #获取左右上下边界

        dLeft=min(left)

        dRight=max(right)

        dTop=max(top)

        dBottom=min(bottom)

       

        #获取分辨率，影像最小分辨率作为数据集分辨率

        dRatioX = min(ratiox)

        dRatioY = min(ratioy)

       

        #计算影像数据集宽度和高度

        nWidth = int((dRight-dLeft)/dRatioX)

        nHeight = int((dTop-dBottom)/dRatioY)

       

        #重新计算，保证分辨率正确

        dRight=dLeft+dRatioX*nWidth

        dBottom=dTop-dRatioY*nHeight

        L = [nWidth, nHeight, dLeft, dTop, dRight, dBottom]

   

    return L

 

def toDB(server, user, pwd, engType, fileType, path):

    files=[]

    files=walkPath(fileType, path)

    print len(files)

    if len(files)>0:

        L=[]

        L = calcDatasetInfo(fileType, files)

        pixType = smu.GetImagePixelFormatName(fileType, files[0])

 

        odsAlias='test'

        if len(L)==6:

            nWidth=L[0]

            nHeight=L[1]

            dLeft=L[2]

            dTop=L[3]

            dRight=L[4]

            dBottom=L[5]

            dtName='test'

            isOpen=smu.OpenDataSource(server,user,pwd, engType, odsAlias)

            smu.DeleteDataset(odsAlias, dtName)

            bCreate = smu.CreateDatasetRaster(odsAlias,dtName,

                    'Image', 'encDCT', pixType,nWidth,nHeight,

                    dLeft, dTop,dRight,dBottom,256)

            writeLog("log.log","calcDatasetInfo 成功")       

            for file in files:

                writeLog("log.log",file+"开始处理")   

                smu.AppendRasterFile(odsAlias,dtName,fileType, file)

                writeLog("log.log",file+"处理完毕\n\n")      

            bBuild=smu.BuildPyramid(odsAlias,dtName)#创建影像金字塔

            if bBuild == 1:

                print "创建金字塔成功"

            else:

                print "创建影像金字塔失败！"

            smu.CloseDataSource(odsAlias)

 

 

#=====================================

def writeLog(logPath, tmpstr):

    time_str = time.strftime("%Y-%m-%d %H:%M:%S ",time.localtime())

    logstr = str(tmpstr) + time_str +'\n'

    print(logstr)

    f = open(logPath, "a")

    f.write(logstr)

    f.close()

 

 

help =u"----------------------------------------------------------\n\

说明:可导入udb或oracle引擎\n\

导入到UDB用法: AppendRasterFile.py ugoPath tif c:/data\n\

导入到Oracle用法: AppendRasterFile.py ugoPath server user pwd tif c:/data\n\

----------------------------------------------------------\n"

 

if __name__=='__main__':

    if len(sys.argv)>2:

        ugo=sys.argv[1]

        if os.path.exists(ugo):

            sys.path.append(ugo)

            import smu

        else:

            print u'组件路径不存在.'

            sys.exit()

    else:

        print help

        sys.exit()

 

    if len(sys.argv) == 4:

        engType='sceUDB'

        fileType=sys.argv[2]

        fileType=getType(fileType)

        path=sys.argv[3]

       

        udb = path+'/test.udb'

        udd = path+'/test.udd'

        if os.path.exists(udb):

            os.remove(udb)

        if os.path.exists(udd):

            os.remove(udd)

 

        print 'toDB ing'

        toDB(udb, '', '', engType, fileType, path)

        smu.Exit()#清空环境，释放内存

    elif len(sys.argv) == 7:

        engType='sceOraclePlus'

        server=sys.argv[2]

        user=sys.argv[3]

        pwd=sys.argv[4]

        fileType=sys.argv[5]

        fileType=getType(fileType)

        path=sys.argv[6]

        toDB(server, user, pwd, engType, fileType, path)

        smu.Exit()#清空环境，释放内存

 

 

 

运行脚本示例：

c:\python27\python.exe d:\p.py "E:\Develop\deskpro6R\Bin" tif d:\tifdata

 

 

Øc:\python27\python.exe版本一定是2.7.3

Ød:\d.py就是拼接GeoTiff或Img格式的SuperMapPy脚本；

Ø"E:\Develop\deskpro6R\Bin"是deskpro安装目录bin，注意确保目录下存在smu.pyd文件；

Øtif指，读取所有GeoTiff；

Ød:\tiffdata是GeoTiff的目录。注意tif文件的后缀名要是小写, 例如d:\tiffdata\1.tif。 

参考：

参考Python for SuperMap  1.0帮助。

来自为知笔记(Wiz)

附件列表