文章详情

短信预约-IT技能 免费直播动态提醒

请输入下面的图形验证码

提交验证

短信预约提醒成功

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

2023-06-02 03:44

关注

  项目需求

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

  项目构想

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

  1、对影像做差值

  2、获取信息发布服务

  项目实现

  观察geoserver文件系统。

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

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

  

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

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

  

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

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

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

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

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

  故利用python 生成四位数不同的随机ID即可。(此时不知道geoserver提供了rest的API)。

  实现流程:

  

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

  代码实现

  项目代码大致分为三个部分(

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

阅读原文内容投诉

免责声明:

① 本站未注明“稿件来源”的信息均来自网络整理。其文字、图片和音视频稿件的所属权归原作者所有。本站收集整理出于非商业性的教育和科研之目的,并不意味着本站赞同其观点或证实其内容的真实性。仅作为临时的测试数据,供内部测试之用。本站并未授权任何人以任何方式主动获取本站任何信息。

② 本站未注明“稿件来源”的临时测试数据将在测试完成后最终做删除处理。有问题或投稿请发送至: 邮箱/279061341@qq.com QQ/279061341

软考中级精品资料免费领

  • 历年真题答案解析
  • 备考技巧名师总结
  • 高频考点精准押题
  • 2024年上半年信息系统项目管理师第二批次真题及答案解析(完整版)

    难度     813人已做
    查看
  • 【考后总结】2024年5月26日信息系统项目管理师第2批次考情分析

    难度     354人已做
    查看
  • 【考后总结】2024年5月25日信息系统项目管理师第1批次考情分析

    难度     318人已做
    查看
  • 2024年上半年软考高项第一、二批次真题考点汇总(完整版)

    难度     435人已做
    查看
  • 2024年上半年系统架构设计师考试综合知识真题

    难度     224人已做
    查看

相关文章

发现更多好内容

猜你喜欢

AI推送时光机
位置:首页-资讯-后端开发
咦!没有更多了?去看看其它编程学习网 内容吧
首页课程
资料下载
问答资讯