利用GDAL實(shí)現(xiàn)柵格影像差值計(jì)算及Geoserver自動(dòng)發(fā)布柵格影像

　　項(xiàng)目需求

　　項(xiàng)目需要對(duì)兩幅柵格影像做差值處理，然后在geoserver上自動(dòng)發(fā)布服務(wù)。

　　項(xiàng)目構(gòu)想

　　仔細(xì)查閱了很多文獻(xiàn)。geoserver上沒有提供直接對(duì)兩幅柵格影像做差值的處理。所以我將步驟分為兩步：

　　1、對(duì)影像做差值

　　2、獲取信息發(fā)布服務(wù)

　　項(xiàng)目實(shí)現(xiàn)

　　觀察geoserver文件系統(tǒng)。

　　就柵格影像的workspace以及styles文件系統(tǒng)為例?！　　　?/p>

　　對(duì)某一指定工作區(qū)添加數(shù)據(jù)存儲(chǔ)以及添加切片。對(duì)某一指定工作區(qū)其namaspace.xml、workspace.xml是相同的。

　　也即是兩個(gè)xml中的namespace與workspace的id生成后不會(huì)因?yàn)樘砑哟鎯?chǔ)而變動(dòng)。

　　對(duì)于柵格影像存儲(chǔ)涉及到coveragestore.xml、layer.xml、coverage.xml　　　　

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

　　從中觀察到各個(gè)文件ID必須不同。涉及到：　　

　　故利用python 生成四位數(shù)不同的隨機(jī)ID即可。(此時(shí)不知道geoserver提供了rest的API)。

　　實(shí)現(xiàn)流程：

　　代碼實(shí)現(xiàn)

　　項(xiàng)目代碼大致分為三個(gè)部分(

　　1、使用gdal對(duì)柵格影像做差值處理

　　2、獲取信息，在相應(yīng)文件系統(tǒng)中自動(dòng)生成xml(分為python生成和geoserver 中rest服務(wù)的api)

　　3、重新啟動(dòng)startup.bat

　　##對(duì)影像處理

　　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 圖像長(zhǎng)度

　　self.bands = self.dataset.RasterCount # TODO 圖像波段數(shù)量

　　self.proj = self.dataset.GetProjection() # todo 地圖投影信息

　　self.geotrans = self.dataset.GetGeoTransform() # todo 仿射矩陣

　　def getRasterInformation(self, nband):

　　band = self.dataset.GetRasterBand(nband) # 獲取波段對(duì)象

　　# 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 新建數(shù)據(jù)集

　　writeable.SetGeoTransform(self.geotrans) # 寫入仿射變換參數(shù)

　　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 給各波段設(shè)置nodata值

　　writeable.GetRasterBand(i + 1).FlushCache() # todo 波段統(tǒng)計(jì)量

　　print(writeable.GetRasterBand(i + 1).GetStatistics(0, 1)) # todo 計(jì)算波段統(tǒng)計(jì)量 輸出為min\max \Mean\stddev

　　def showImage(self, r, g, b):

　　img2 = cv2.merge([r, g, b])

　　plt.imshow(img2)

　　plt.xticks([]), plt.yticks([]) # 不顯示坐標(biāo)軸

　　plt.title("CHA")

　　plt.show()

　　## 生成xml 以layer.xml為例

　　import xml.dom.minidom

　　import random

　　import os 　無錫人流醫(yī)院哪家好 http://www.wxbhnkyy120.com/

　　def writelayerXml(self):

　　fp = open(self.mainpath + self.vapath + self.vapath + "\layer.xml", 'w')

　　# TODO 在內(nèi)存中創(chuàng)建一個(gè)空的文檔

　　doc = xml.dom.minidom.Document()

　　# TODO 創(chuàng)建一個(gè)根節(jié)點(diǎn)Managers對(duì)象

　　root = doc.createElement('layer')

　　# todo 將根節(jié)點(diǎn)添加到文檔對(duì)象中

　　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")

　　###啟動(dòng)服務(wù)

　　import win32api

　　def startbat():

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

　　執(zhí)行完成之后，即可看到geoserver中新增了數(shù)據(jù)存儲(chǔ)以及服務(wù)。

　　其中關(guān)于geoserver rest服務(wù)中的API下一節(jié)再講述。