利用GDAL实现栅格影像差值计算及Geoserver自动发布栅格影像

项目需求

项目需要对两幅栅格影像做差值处理，然后在geoserver上自动发布服务。

项目构想

仔细查阅了很多文献。geoserver上没有提供直接对两幅栅格影像做差值的处理。所以我将步骤分为两步：
1、对影像做差值
2、获取信息发布服务
关于自动发布的思想借鉴于这位博主的思想，可查阅这篇博客
https://www.cnblogs.com/naaoveGIS/p/4212093.html

项目实现

观察geoserver文件系统。
就栅格影像的workspace以及styles文件系统为例。

对某一指定工作区添加数据存储以及添加切片。对某一指定工作区其namaspace.xml、workspace.xml是相同的。

也即是两个xml中的namespace与workspace的id生成后不会因为添加存储而变动。

对于栅格影像存储涉及到coveragestore.xml、layer.xml、coverage.xml


其中coverage篇幅较大，较不同地方在于

从中观察到各个文件ID必须不同。涉及到：

故利用python 生成四位数不同的随机ID即可。（此时不知道geoserver提供了rest的API）。
实现流程：

代码实现

项目代码大致分为三个部分(
1、使用gdal对栅格影像做差值处理
2、获取信息，在相应文件系统中自动生成xml（分为python生成和geoserver 中rest服务的api)
3、重新启动startup.bat

##对影像处理
import gdal, gdalconst, numpy
import cv2
import matplotlib.pyplot as plt
class ReadRaster:
    def __init__(self, path, ):
        self.dataset = gdal.Open(path, gdal.GA_ReadOnly)
        self.rows = self.dataset.RasterYSize  # todo  图像宽度
        self.cols = self.dataset.RasterXSize  # todo  图像长度
        self.bands = self.dataset.RasterCount  # TODO 图像波段数量
        self.proj = self.dataset.GetProjection()  # todo 地图投影信息
        self.geotrans = self.dataset.GetGeoTransform()  # todo 仿射矩阵

    def getRasterInformation(self, nband):
        band = self.dataset.GetRasterBand(nband)  # 获取波段对象
        # data = band.ReadAsArray(0, 0, self.cols, self.rows).astype(numpy.float)  #获取波段信息
        data = band.ReadAsArray(0, 0, self.cols, self.rows)  # 获取波段信息
        return data

    def writeRasterInformation(self, data, Savepath, nband):
        driver = self.dataset.GetDriver()
        writeable = driver.Create(Savepath, self.cols, self.rows, self.bands, gdal.GDT_Byte)  # TODO  新建数据集
        writeable.SetGeoTransform(self.geotrans)  # 写入仿射变换参数
        writeable.SetProjection(self.proj)  # 写入投影
        for i in range(nband):
            writeable.GetRasterBand(i + 1).WriteArray(data[i], 0, 0)
            writeable.GetRasterBand(i + 1).SetNoDataValue(0)  # todo 给各波段设置nodata值
            writeable.GetRasterBand(i + 1).FlushCache()  # todo 波段统计量
            print(writeable.GetRasterBand(i + 1).GetStatistics(0, 1))  # todo 计算波段统计量  输出为min\max \Mean\stddev

    def showImage(self, r, g, b):
        img2 = cv2.merge([r, g, b])
        plt.imshow(img2)
        plt.xticks([]), plt.yticks([])  # 不显示坐标轴
        plt.title("CHA")
        plt.show()


## 生成xml 以layer.xml为例

import xml.dom.minidom
import random
import os

    def writelayerXml(self):
        fp = open(self.mainpath + self.vapath + self.vapath + "\layer.xml", 'w')
        # TODO  在内存中创建一个空的文档
        doc = xml.dom.minidom.Document()
        # TODO 创建一个根节点Managers对象
        root = doc.createElement('layer')
        # todo 将根节点添加到文档对象中
        doc.appendChild(root)
        # id\name\description\type\enabled\workspace\__default\url
        nodeid = doc.createElement('id')
        nodename = doc.createElement('name')
        nodetype = doc.createElement('type')
        nodedefaultStyle = doc.createElement('defaultStyle')
        nodedefaultStyleid = doc.createElement('id')
        noderesource = doc.createElement('resource')
        noderesource.setAttribute('class', 'coverage')
        noderesourceid = doc.createElement('id')
        nodeattribution = doc.createElement('attribution')
        nodelogoWidth = doc.createElement('logoWidth')
        nodelogoHeight = doc.createElement('logoHeight')

        nodename.appendChild(doc.createTextNode(self.rastername))
        root.appendChild(nodename)
        nodeid.appendChild(doc.createTextNode(ProduceXml.LayerInfoImpl+self.password3))
        root.appendChild(nodeid)
        nodetype.appendChild(doc.createTextNode('RASTER'))
        root.appendChild(nodetype)
        nodedefaultStyleid.appendChild(doc.createTextNode(ProduceXml.StyleInfoImpl))
        nodedefaultStyle.appendChild(nodedefaultStyleid)
        root.appendChild(nodedefaultStyle)
        noderesourceid.appendChild(doc.createTextNode(ProduceXml.CoverageInfoImpl+self.password2))
        noderesource.appendChild(noderesourceid)
        root.appendChild(noderesource)
        nodelogoWidth.appendChild(doc.createTextNode('0'))
        nodelogoHeight.appendChild(doc.createTextNode('0'))
        nodeattribution.appendChild(nodelogoWidth)
        nodeattribution.appendChild(nodelogoHeight)
        root.appendChild(nodeattribution)

        doc.writexml(fp, indent='\t', addindent='\t', newl='\n', encoding="utf-8")



###启动服务

import win32api
def startbat():
    win32api.ShellExecute(0, 'open', r'C:\Program Files (x86)\GeoServer 2.14.2\bin\startup.bat', '', '', 1)

执行完成之后，即可看到geoserver中新增了数据存储以及服务。

其中关于geoserver rest服务中的API下一节再讲述。