python 经纬度 内切圆

Python GDAL灰度图像转RGB图像

工作环境：Python3.6 GDAL Numpy
工作时间：2019/03/10-2019/03/17
本文在AmosHawk,WHU,LIESMARS的代码的基础上进行改进。源代码支持将图像按指定长度切割，但是切割之后的图像不能按照切割位置写入切割后的地理坐标。

文章目录

• Python GDAL灰度图像转RGB图像

• 1.处理思路
• 2.代码结构
• 3.每个模块分析

1.处理思路

1.判断任务类型：输入为文件夹还是单个文件，将多个文件按单个图像循环交给处理函数处理；
1.处理函数首先读取原图像xy像素数量，和目像素数量对比，确定切割数量；
2.按要求依次取出每个目标图片的各个通道的数值；
3.将目标图像的像素的值和图像的基本信息依次写给输出图像。

2.代码结构

与灰度图像转RGB图像一样，由几个def函数和一个类及其方法构成。

import...
def OpenArray( array, prototype_ds = None,xoff=0, yoff=0):...
def resizesingefile(singlefileName, outputFile, heightsubImage, widthsubImage):
class cut_img:
    def resize(self, mode, imgtype, inputfile, outputfile, heightsubImage, widthsubImage):

3.每个模块分析

1.OpenArray
gdal_array.OpenArray(array, prototype_ds=None)结合prototype_ds的图像信息（如地理坐标等）和array的图像像素（如RGB的像素值）数组并赋给返回对象。但是不能改写返回图像的地理坐标。（事实上，切割图片的输出和输入图像信息只有地理坐标需要修改）因此，改写该函数，命名为OpenArray，该函数多了两个输入值，可以将该值指定返回对象的地理坐标。

def OpenArray( array, prototype_ds = None,xoff=0, yoff=0):
    ds = gdal_array.OpenNumPyArray( array )
    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdal_array.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff)
    return ds

2.resizesingefile单个图像切割
类的方法会调用该函数，该函数接受的文件是无差别的所有格式的文件。因此需要函数可以辨别哪些文件可以被处理。这里我用了一个try…except…+while True…break…的方法。gdalnumeric.LoadFile如果无法加载文件，说明gdal不支持该文件类型，也就会报错，从而执行break接触函数。
图像分割思路：
我们假设有一个图片会被切割成2x2的4张小图像，则依次从原图中取出(1,1),(1,2),(2,1),(2,2)的每个像素点上的n通道的值（例如3bandsRGB图像）组成一个(x,y,3)的三维数组，将这四个三维数组依次写入一个数组subImages中，组成一个四维数组。
从这个四位数组中依次将每个(x,y,3)的像素和源文件的图像信息，外带xy地理坐标通过上面的OpenArray函数合并。

def resizesingefile(singlefileName, outputFile, heightsubImage, widthsubImage):
    while True:
        # Load the source data as a gdalnumeric array
        try:
            srcArray = gdalnumeric.LoadFile(singlefileName)
        except:
            break
        # Also load as a gdal image to get geotransform
        # (world file) info
        srcImage = gdal.Open(singlefileName)
        if srcImage is None:
            print('can not open', singlefileName, 'with gdal')
        oriHei = srcImage.RasterYSize
        oriWid = srcImage.RasterXSize
        oriBandNum = srcImage.RasterCount
        oriTransform = srcImage.GetGeoTransform
        if (oriHei <= heightsubImage | oriWid <= widthsubImage):
            print('image :', singlefileName, 'is smaller than the specified shape')
        # creat a new subcontent to store the subimages and place it to the upper content
        # 在输入目录下创建目标文件夹
        dirname, filename = os.path.split(singlefileName)
        filename = filename.split('.')
        newSubContent = outputFile + '//' + filename[0]
        if (os.path.exists(newSubContent) == False):
            os.mkdir(newSubContent)       
        # calculate the numbers by row and coloum by the specific width and heigh
        nRowNums = ceil(oriHei / heightsubImage)
        nColNums = ceil(oriWid / widthsubImage)
        # build a list to store the subimage data for the moment
        subImages = []
        x = []
        y = []
        # begin to crop the image
        # 分单通道和多通道
        if oriBandNum == 1:
            for i in range(0, nRowNums):
                for j in range(0, nColNums):
                    # subImage = oriImage[i*heightsubImage:(i+1)*heightsubImage,j*widthsubImage:(j+1)*widthsubImage]
                    clip = srcArray[i * heightsubImage:(i + 1) * heightsubImage, j * widthsubImage:(j + 1) * widthsubImage]
                    subImages.append(clip)
                    # 经纬度
                    y.append(i * heightsubImage)
                    x.append(j * widthsubImage)
        else:
            for i in range(0, nRowNums):
                for j in range(0, nColNums):
                    # subImage = oriImage[i*heightsubImage:(i+1)*heightsubImage,j*widthsubImage:(j+1)*widthsubImage]
                    clip = srcArray[:, i * heightsubImage:(i + 1) * heightsubImage,
                           j * widthsubImage:(j + 1) * widthsubImage]
                    subImages.append(clip)
                    y.append(i * heightsubImage)
                    x.append(j * widthsubImage)
        # wirte the image to the new created subcontent
        for j in range(1, len(subImages) + 1):
            # print('begin to write :', j, 'th subimage of', file)
            savefile = newSubContent + "//" + filename[0] + '_' + str(j) + '.' + filename[1]
            gtiffDriver = gdal.GetDriverByName('GTiff')
            if gtiffDriver is None:
                raise ValueError("Can't find GeoTiff Driver")
            gtiffDriver.CreateCopy(savefile,
                                   OpenArray(subImages[j - 1], prototype_ds=singlefileName, xoff=x[j - 1],
                                             yoff=y[j - 1]))
            print('finish writting:',savefile)
        break

3.类：接受输入并分析任务，将任务分解成单个文件任务交给函数resizesingefile执行

class cut_img:
    def resize(self, mode, imgtype, inputfile, outputfile, heightsubImage, widthsubImage):
        time_start = time.time()
        if mode == [1]:
            resizesingefile(inputfile, outputfile, heightsubImage, widthsubImage)
        elif mode == [2]:
            for file in os.listdir(inputfile):
                singlefileName = inputfile+"\\"+file
                singlefileForm = os.path.splitext(singlefileName)[1][1:]
                if(singlefileForm == imgtype):
                    resizesingefile(singlefileName, outputfile, heightsubImage, widthsubImage)
                elif (imgtype == 'all'):
                    resizesingefile(singlefileName, outputfile, heightsubImage, widthsubImage)
        time_end = time.time()
        return time_end - time_start

import如下：

import os
from math import ceil
from osgeo import gdal, gdalnumeric,gdal_array
import time