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

项目需求

创新互联建站-专业网站定制、快速模板网站建设、高性价比莘县网站开发、企业建站全套包干低至880元,成熟完善的模板库,直接使用。一站式莘县网站制作公司更省心,省钱,快速模板网站建设找我们，业务覆盖莘县地区。费用合理售后完善，10余年实体公司更值得信赖。

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

项目构想

仔细查阅了很多文献。geoserver上没有提供直接对两幅栅格影像做差值的处理。所以我将步骤分为两步：

1、对影像做差值

2、获取信息发布服务

项目实现

观察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 无锡人流医院哪家好 http://www.wxbhnkyy120.com/

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下一节再讲述。