森林火险气象预报图GIS应用

背景

中央气象台每日发布一幅全国范围内的森林火险气象预报图，基于森林火险气象等级分为三类：较高、高、极高。并且对相关区域进行了文字说明。在公司的自然灾害监测预警平台中，可以将这部分官方数据加载进去，以达到针对火险预警的指示性作用，结合现有的火点实时动态显示，进一步完善火情监测预报的功能体系。

在上述功能需求下，汇总了以下的问题和解决方案。数据获取、数据处理、数据展示。数据来源于网站，需要进行爬取；图片没有坐标系统，需要进行配准；只需要火险等级区域，需要进行信息提取；提取后数据是单波段栅格数据，需要进行颜色匹配；前端动态加载，需要使用动态服务。

数据获取
1. 数据来源

中央气象台森林火险气象预报：http://www.weather.com.cn/index/zxqxgg1/slqxyb.shtml

1. 数据内容

数据当天下午6点左右发布第二天的预报内容，包括标题、单位、发布时间、预报范围信息及图片。

1. 数据爬取
  1. 流程

1. 1. Python脚本

from bs4 import BeautifulSoup

from lxml import html

import openpyxl

from openpyxl.styles import PatternFill, Font, Alignment, Border, Side

import requests

import urllib.request

import re

import datetime

import time

import threading

import os

import cx_Oracle

import pandas as pd

from sqlalchemy import create_engine

import numpy as np

excelName = '森林火险预报记录.xlsx'

ExcelFullName = os.path.join(os.getcwd(), excelName)

os.environ['NLS_LANG'] = 'AMERICAN_AMERICA.AL32UTF8' #保证select出来的中文显示没有问题

def createExcel():

if not os.path.exists(ExcelFullName):

wbExport = openpyxl.Workbook()

wbExport.active

ws1 = wbExport["Sheet"]

ws1.title = '记录' # 重命名为站点名

tableTitle = ['Date', 'Title', 'Time', 'Content', 'Picture', 'Remarks']

for col in range(len(tableTitle)):

c = col + 1

ws1.cell(row=1, column=c).value = tableTitle[col]

wbExport.save(filename=ExcelFullName)

def saveSLQX():

createExcel()

wb = openpyxl.load_workbook(ExcelFullName)

ws0 = wb["记录"]

inputValues = ['', '', '', '', '', '']

today = datetime.date.today()

formatted_today = today.strftime('%Y%m%d')

imgname = os.path.join(

os.getcwd(), 'sevp_nmc_'+formatted_today+'.jpg')

url = "http://www.weather.com.cn/index/zxqxgg1/slqxyb.shtml"

page = urllib.request.urlopen(url).read().decode("utf-8", "ignore")

divtitle = re.compile(r'<div class="title">(.*?)</div>')

divpageContent = re.compile(r'<div class="pageContent"><p style="text-indent:2em;">(.*?)</p></div>')

imgRe = re.compile(r'src="(.*?\.jpg)"')

titles = divtitle.findall(page)

contents = divpageContent.findall(page)

images = imgRe.findall(page)

imgs = []

for img in images:

imgs.append(img)

comp = re.compile(r'<[^>]+>', re.S)

pTitle = comp.sub(' ', titles[0]).strip()

pContent = comp.sub('_', contents[0])

pTitle1 = pTitle.split(' ')

inputValues[0] = formatted_today

inputValues[1] = pTitle1[0]

inputValues[2] = pTitle1[2]+' '+pTitle1[3]

pContent1 = pContent.split('_')[0]

pContent1 = pContent1+'。'

pContent2 = pContent.split('_')[1:-1]

pRemark = pContent.split('_')[-1]

for c in pContent2:

pContent1 = pContent1+c

inputValues[3] = pContent1.replace(' ', '')

inputValues[4] = imgname

inputValues[5] = pRemark

# inputValues.append(pRemark)

ws0.append(inputValues)

# 列宽

for column in ws0.columns:

if(column[00].column_letter == 'A'):

ws0.column_dimensions[column[00].column_letter].width = 10

elif(column[00].column_letter == 'B'):

ws0.column_dimensions[column[00].column_letter].width = 30

elif(column[00].column_letter == 'C'):

ws0.column_dimensions[column[00].column_letter].width = 20

elif(column[00].column_letter == 'D'):

ws0.column_dimensions[column[00].column_letter].width = 50

elif(column[00].column_letter == 'E'):

ws0.column_dimensions[column[00].column_letter].width = 30

elif(column[00].column_letter == 'F'):

ws0.column_dimensions[column[00].column_letter].width = 30

else:

pass

# 行高

for row in ws0.rows:

for cell in row:

cell.alignment = Alignment(vertical='center') # 垂直居中

row[3].alignment = Alignment(vertical='center', wrapText=True)

row[4].alignment = Alignment(vertical='center', wrapText=True)

row[5].alignment = Alignment(vertical='center', wrapText=True)

response = urllib.request.urlretrieve(url=imgs[0], filename=imgname)

wb.save(filename=ExcelFullName)

now_time = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")

insert2Table(inputValues)

print('成功！'+now_time)

def insert2Table(listvalues):

wholeValues=""

for lv in listvalues:

wholeValues=wholeValues+",'"+lv+"'"

values=wholeValues[1:]

# conn = cx_Oracle.connect('sde/sde@122.112.200.75/orcl')

conn = cx_Oracle.connect('SLHQ/SLHQ@localhost/orcl')

cursor = conn.cursor()

query = "INSERT INTO HXDJ(dateymd,title,time,content,picture,remarks) VALUES ({})"

query=query.format(values)

cursor.execute(query)

conn.commit()

cursor.close()

conn.close()

result = saveSLQX()

1. 1. 代码详解

代码基于python3.X进行编写。引入的函数库有BeautifulSoup4，lxml，openpyxl，requests，urllib，cx_Oracle，pandas，sqlalchemy，numpy。

createExcel()函数用于创建excel文件，用于之后爬取数据的文字内容存储。

saveSLQX()函数用于图片下载另存，文字内容编辑插入到excel表格。

insert2Table(listvalues)函数基于文字内容，将内容插入到Oracle数据库对应的表格中。

由于Oracle数据库插入数据时，id没有办法自增，因此，在创建的表格后，还需创建序列和触发器。这样在插入语句中就不用增加id字段，从而实现数据有序存储到Oracle中。

数据处理
1. 处理工具

ArcGIS提供一套自身的py库，arcpy。在对图片、影像等处理中，arcpy有着专业性的优势。此处用到的工具有影像配准，栅格数据提取，重分类，添加色彩映射表。

1. 流程

Arcmap中使用model builder进行流程化操作，然后将其转换成python文件，再对python进行修改，可以快捷、灵活、准确的完成脚本的编辑。

1. 问题

在过年及疫情期间，脚本一直正常运行，下载数据，提取数据，但是有些图像提取的结果为空，查看原图片，是有等级范围的。在比较原始图片时，发现其数据的物理存储大小不一样，但是其分辨率大小一致。通过原图加载发现，一张图片颜色的阈值是16种，另一张图片颜色的阈值是256种。如下图所示。

出现上述原因应该是原图的输出模式上，其色彩映射表的选择上是4bit unsigned int(0~15)与8bit unsigned int (0~255)的原因。因此，在对数据提取时，对应的代码进行判断，并且增加一个重分类的步骤。因为在16bit颜色值下，其高等级范围内含两种颜色值。如下图所示。

Figure: 8bit颜色值

Figure: 4bit颜色值

根据发现的问题，结合图像最终的展示效果需求，在图片重分类后，加上4bit原图的那种色彩映射表后输出即可生成有颜色的tif图片，便于直接使用。

1. Python脚本

import arcpy

import datetime

import os

try:

osstart = "G:\DATAALL\HxdjImages"

arcpy.env.workspace = osstart

arcpy.env.overwriteOutput = True

print(osstart)

formatted_today = datetime.date.today().strftime('%Y%m%d')

YuanImgname = os.path.join(osstart, 'sevp_nmc_'+formatted_today+'.jpg')

if not os.path.exists(YuanImgname):

print(r'Image: "'+YuanImgname+'" is not exist!')

quit()

OutImgFolder = os.path.abspath(os.path.join(osstart, 'outmap'))

OutImgfile = os.path.join(

OutImgFolder, 'outmap_sevp_nmc_'+formatted_today+'.tif')

print(OutImgfile)

size = os.path.getsize(YuanImgname)/1024

# Check out any necessary licenses

arcpy.CheckOutExtension("spatial")

# Script arguments

inname = ''

if inname == '#' or not inname:

inname = YuanImgname # provide a default value if unspecified

resul_inamet_tif = ''

if resul_inamet_tif == '#' or not resul_inamet_tif:

resul_inamet_tif = OutImgfile # provide a default value if unspecified

# Local variables:

geomap = arcpy.CreateScratchName(

"temp", data_type="RasterDataset", workspace=arcpy.env.workspace)

outmap_inname_tif=arcpy.CreateScratchName(

"tempr", data_type="RasterDataset", workspace=arcpy.env.workspace)

peizhun_txt = r"G:\DATAALL\HxdjImages\tool\project.txt"

outmap_sevp_nmc_hxdj_clr = r"G:\\DATAALL\\HxdjImages\\tool\\outmap_sevp_nmc_hxdj.clr"

# Process: Warp From File

tempEnvironment0 = arcpy.env.outputCoordinateSystem

arcpy.env.outputCoordinateSystem = "PROJCS['Asia_Lambert_Conformal_Conic',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',105.0],PARAMETER['Standard_Parallel_1',30.0],PARAMETER['Standard_Parallel_2',62.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"

arcpy.WarpFromFile_management(

inname, geomap, peizhun_txt, 'SPLINE', 'NEAREST')

arcpy.env.outputCoordinateSystem = tempEnvironment0

# Process: Extract by Attributes

tempEnvironment0 = arcpy.env.mask

arcpy.env.mask = r"G:\DATAALL\HxdjImages\tool\mask.shp"

condi = "\"Value\" =180 OR \"Value\" =198 OR \"Value\" =204 OR \"Value\" =210"

if size < 400:

condi = "\"Value\" =8 OR \"Value\" =9 OR \"Value\" =10"

arcpy.gp.ExtractByAttributes_sa(

geomap, condi, outmap_inname_tif)

arcpy.env.mask = tempEnvironment0

# Process: Reclassify

arcpy.gp.Reclassify_sa(outmap_inname_tif, "VALUE", "180 8;198 9;204 9;210 10", resul_inamet_tif, "DATA")

# Process: Add Colormap

arcpy.AddColormap_management(resul_inamet_tif, "", outmap_sevp_nmc_hxdj_clr)

arcpy.Delete_management(geomap)

arcpy.Delete_management(outmap_inname_tif)

print("Job completed!")

except:

print(arcpy.GetMessages())

1. 代码详解

代码是基于arcpy进行编写的，由于arcmap是32位的软件，其自动配置的是python2.7，语法与python2.X相同。

其中，涉及到数据操作处理，均是arcpy调用相应函数，给定参数即可初步实现，在处理过程中，会增加中间变量，在数据处理完成后，删除即可。

脚本定时自动执行
1. 任务计划程序库

计算机管理→系统工具→任务计划程序→任务计划程序库→新建任务

根据要求进行任务创建，任务常规、触发器、操作对象、条件、相关设置等。

1. 案例

以森林火险气象预报图及内容的爬取脚本为例。

1、常规设置

2、触发器，设置定时器

3、设置启动程序

启动程序位置，执行的脚本路径如上图所示。

4、其他（启动条件、设置）

注意：由arcpy编写的脚本，其需要使用arcmap安装路径下的python打开，并且由于python2.X对中文路径不友好，需要注意脚本中的路径变量。

数据展示

生成的图片存放在服务器中，对于每天都生成的图片，每天发布成服务是不合理的。arcgis server提供动态服务功能，即可以将shapefile文件夹、本地数据库、空间数据库、raster文件夹注册到服务器，支持动态的调用注册文件夹内的数据。

Figure: 注册动态数据空间

function initTif() {

var d = new Date();

var curr_date = d.getDate()-1;

var curr_month = d.getMonth() + 1;

var curr_year = d.getFullYear();

String(curr_month).length < 2 ? (curr_month = "0" + curr_month) : curr_month;

String(curr_date).length < 2 ? (curr_date = "0" + curr_date) : curr_date;

var yyyyMMdd = curr_year + "" + curr_month + "" + curr_date;

tifName = "outmap_sevp_nmc_" + yyyyMMdd;

addHanzaiTif(tifName);

}

function addHanzaiTif(result) {

map.remove(gh1);

gh1 = null;

gh1 = new MapImageLayer({

url: "http://122.112.200.75:6080/arcgis/rest/services/hanzai/HXDJServer/MapServer",

sublayers: [

{

id: [0],

opacity: 0.8,

visible: true,

source: {

type: "data-layer",

dataSource: {

type: "raster",

workspaceId: "hxdjRaster",//outraster

dataSourceName: result + ".tif"

}

});

map.add(gh1);//加载图层

}

上述代码为调用当天的森林火险气象预报数据，其数据展示效果如下图所示。

森林火险气象预报图GIS应用

天立泰在线知识库

标签云

友情链接

网站信息