付嬈

（遼寧省自然資源事務(wù)服務(wù)中心，遼寧 沈陽 110032）

焚燒秸稈會(huì)在短時(shí)間內(nèi)嚴(yán)重影響到空氣質(zhì)量，污染城市大氣環(huán)境。2017年10月9日，哈爾濱市人民政府印發(fā)《哈爾濱市禁止野外焚燒秸稈改善大氣環(huán)境質(zhì)量實(shí)施方案》，根據(jù)《方案》，哈爾濱將建立市級(jí)領(lǐng)導(dǎo)包區(qū)縣(市)、區(qū)縣(市)干部包鄉(xiāng)(鎮(zhèn))、鄉(xiāng)(鎮(zhèn))干部包村、鄉(xiāng)(鎮(zhèn))村(屯)干部包農(nóng)戶和地塊的逐級(jí)包保責(zé)任制，形成以鄉(xiāng)(鎮(zhèn))、村(屯)、農(nóng)戶、地塊為單元的網(wǎng)格化管理體系。快速獲取火點(diǎn)數(shù)據(jù)并對(duì)秸稈焚燒發(fā)生的地點(diǎn)進(jìn)行實(shí)時(shí)核查，成為了落實(shí)逐級(jí)包保責(zé)任制關(guān)鍵的一環(huán)。

一、研究區(qū)和數(shù)據(jù)來源

哈爾濱市位于東北平原東北部地區(qū)、黑龍江省南部，地勢(shì)東南高,西北低,根據(jù)2016年哈爾濱市統(tǒng)計(jì)年鑒，2016年年末哈爾濱市國土面積5.31萬平方公里，耕地面積1.97756萬平方公里。約占總面積的37.24%。其中，市區(qū)耕地面積為0.540434萬平方公里，市轄縣（市）耕地面積為1.437131萬平方公里。

所用的火點(diǎn)經(jīng)緯度數(shù)據(jù)來自網(wǎng)站https://firms.modaps.eosdis.nasa.gov/，該網(wǎng)站一共提供了3種不同的數(shù)據(jù)格式，包括Google Earth KML，CSV表格，Shapefile；并且每種數(shù)據(jù)格式都提供了3種不同時(shí)間段的數(shù)據(jù)，包括24小時(shí)內(nèi)的火點(diǎn)數(shù)據(jù)，48小時(shí)的火點(diǎn)數(shù)據(jù)和最近7天內(nèi)的火點(diǎn)數(shù)據(jù)。

其中，KML （Keyhole Markup Language），是一個(gè)基于XML語法和文件格式的文件，用來描述和保存地理信息如點(diǎn)、線、圖片、折線并在Google Earth客戶端顯示；CSV（Comma-Separated Values）即逗號(hào)分隔值，是一種通用的、相對(duì)簡(jiǎn)單的文件格式，其文件以純文本形式存儲(chǔ)表格數(shù)據(jù)（數(shù)字和文本）；Shapefile 文件是描述控件數(shù)據(jù)的幾何和屬性特征的非拓?fù)鋵?shí)體矢量數(shù)據(jù)結(jié)構(gòu)的一種格式，一個(gè)Shapefile文件至少包括三個(gè)文件：一是主文件（*.shp），存儲(chǔ)地理要素的幾何圖形的文件；二是索引文件（*.shx），存儲(chǔ)圖形要素與屬性信息索引的文件；三是dBASE表文件（*.dbf），即存儲(chǔ)要素信息屬性的dBase表文件。

不同的數(shù)據(jù)格式，但是都用相同的字段標(biāo)識(shí)，文件中共包含13個(gè)字段：經(jīng)度(Latitude)，緯度(Longitude)，亮度溫度I-4(Bright_ti4)，沿掃描像素大小(Scan)，沿著軌道像素大小(Track)，采集日期(Acq_Date)，采集時(shí)間(Acq_Time)， 衛(wèi) 星 (Satellite)， 置 信 度 (Confidence)， 版 本(Version)，亮度溫度I-5(Bright_ti5)，火輻射能量(FRP)，白天或晚上(DayNight)。本文以爬取24小時(shí)內(nèi)的Shapefile數(shù)據(jù)為例。

二、Python語言數(shù)據(jù)處理

Python是一門解釋型語言，使用專門的解釋器對(duì)源程序進(jìn)行逐行解釋成某個(gè)特定平臺(tái)的機(jī)器碼，在運(yùn)行時(shí)將程序翻譯成機(jī)器語言，解釋型語言相當(dāng)于把編譯型語言中的編譯和解釋過程混合在一起同時(shí)完成。它也常被稱作膠水語言，能夠把用其他語言（如C，C++等）制作的各種模塊，引入到腳本中使用[1]。Python具有豐富的開源庫，可擴(kuò)展性與可嵌入性強(qiáng)，代碼簡(jiǎn)單易懂、可讀性強(qiáng)等特點(diǎn)。

Python有可定義的第三方庫可以使用，如ArcGIS中自帶的ArcPy庫專門用于地理處理；包括正則表達(dá)式、文檔生成、單元測(cè)試、線程、數(shù)據(jù)庫、網(wǎng)頁瀏覽器、CGI、FTP、電子郵件、XML、XML-RPC、HTML、WAV文件、密碼系統(tǒng)、GUI（圖形用戶界面）、Tk和其他與系統(tǒng)有關(guān)的操作。除了Python標(biāo)準(zhǔn)庫以外，還可以引用許多高質(zhì)量的庫，如wxPython、Twisted和Python圖像庫等等。還可以在程序中引入深度學(xué)習(xí)的各種框架，如TensorFlowKeras，Caffe，MXNet等等。

1.原始數(shù)據(jù)處理

（1）數(shù)據(jù)下載

Python標(biāo)準(zhǔn)庫中的urllib2模塊，定義了一些類和方法主要用于實(shí)現(xiàn)對(duì)HTTP通信協(xié)議的支持，urllib2支持HTTP代理、HTTP簡(jiǎn)單認(rèn)證、跳轉(zhuǎn)、Cookie等功能。urllib2模塊還支持對(duì)HTTP請(qǐng)求報(bào)文的頭和實(shí)體進(jìn)行增改，對(duì)HTTP應(yīng)答報(bào)文的頭和正文進(jìn)行讀取[2]。文中利用urllib2模塊下載Shapefile數(shù)據(jù)格式的壓縮包數(shù)據(jù)，并設(shè)置數(shù)據(jù)存儲(chǔ)的位置。導(dǎo)入urllib2庫下載壓縮包，并導(dǎo)入Python標(biāo)準(zhǔn)庫中的os模塊，用于程序中切換文件路徑、自動(dòng)創(chuàng)建或者移除文件夾、更改文件的名字、判斷指定文件夾在指定路徑下是否存在等操作。

導(dǎo)入Python的標(biāo)準(zhǔn)庫time模塊獲取當(dāng)前日期與時(shí)間，并轉(zhuǎn)換日期的格式。

①獲取系統(tǒng)時(shí)間，并將日期設(shè)為全局變量，以日期命名所有下載的文件，以便以后查詢相應(yīng)時(shí)間段的數(shù)據(jù)：

import urllib2,os,time

global _date

now = int(time.time())

timeArray = time.localtime(now)

_date = time.strftime("%Y%m%d", timeArray)

②new_path = "D:\data"

_pathZipFolder = new_path + "\"+zip+"\"

url = 'https://firms.modaps.eosdis.nasa.gov/data/active_fire/c6/shapes/zips/MODIS_C6_Russia_and_Asia_24h.zip'

# header變量是一個(gè)包含各種瀏覽器信息的字典，目的是用來偽裝瀏覽器爬取數(shù)據(jù)

header = {'User-Agent':'Mozilla/5.0 (Windows NT 6.1;WOW64) AppleWebKit/537.36 (KHTML, like Gecko)Chrome/31.0.1650.63 Safari/537.36'}

os.chdir(_pathZipFolder)

req = urllib2.Request(url,headers = header)

resp = urllib2.urlopen(req)

data = resp.read()

save_name = _date +"_" +"24h"+".zip"

with open(save_name, "wb") as code:

code.write(data)

上述代碼在程序中下載的壓縮包在計(jì)算機(jī)中的排列方式如圖1。

（2）數(shù)據(jù)解壓以及重命名

導(dǎo)入解壓zip數(shù)據(jù)格式的第三方Python庫zipfile將前一步驟中生成的壓縮包解壓到指定文件路徑下，利用for循環(huán)將文件夾中的所有壓縮包解壓并更改解壓后原來文件的命名，以日期加上下載數(shù)據(jù)的時(shí)間段命名的方式代替,方便以后研究使用對(duì)應(yīng)日期的數(shù)據(jù)，如20180301_24h.shp。由于解壓后的數(shù)據(jù)沒有prj格式的數(shù)據(jù)，即Shapefile數(shù)據(jù)的空間參考文件，所以最后解壓完成后需要將研究區(qū)耕地矢量數(shù)據(jù)中的空間參考文件拷貝到解壓后的指定目錄下，并且文件的命名方式與解壓更改的命名保持一致，如20180301_24h.prj，拷貝文件需要用到Python自帶的文件、文件夾壓縮包的處理模塊shutil：

①定義變量作為解壓后文件的存儲(chǔ)路徑：

_pathShpFolder=new_path + "\"+"shp"+"\"

_pathNewDate=_pathShpFolder+"24h"+"\"

②判斷文件路徑是否已經(jīng)存在，若不存在，則創(chuàng)建；若存在，則將解壓后的文件放置在該文件目錄下：

isExists=os.path.exists(_pathNewDate)

if not isExists:

os.makedirs(_pathNewDate)

print _pathNewDate+' directory made successfully'

else:

for file in os.listdir(_pathNewDate):

targetFile = os.path.join(_pathNewDate, file)

if os.path.isfile(targetFile):

os.remove(targetFile)

zipFile=_pathZipFolder+_date+"_"+"24h" +".zip"

③切換當(dāng)前路徑到shp文件存儲(chǔ)的路徑下，并將耕地的空間參考信息prj文件拷貝到路徑下：

os.chdir(_pathNewDate)

shutil.copyfile(_pathPrjFile,"prj.prj")

圖1 壓縮包在計(jì)算機(jī)中的排列方式

newShpDir = os.listdir(_pathNewDate)

④利用for循環(huán)，提取文件夾下所有數(shù)據(jù)的擴(kuò)展名，在不改變文件自身的屬性的條件下，改變文件命名方式：

for temp in newShpDir:

suffix = os.path.splitext(temp)[1]

new_name = _date+"_" +"24h" + suffix

os.rename(temp,new_name)

上述代碼在程序中解壓壓縮包后的文件在計(jì)算機(jī)中的排列方式如圖2。

2.地理數(shù)據(jù)處理

ArcPy的制圖模塊可以根據(jù)ArcMap會(huì)話中的工作流進(jìn)行工作，一個(gè)典型的工作流程包括：打開地圖文檔，修改數(shù)據(jù)框?qū)傩裕虞d圖層，修改圖層屬性，編輯頁面布局中的元素，將地圖文檔導(dǎo)出成圖片格式。上述過程都可在腳本中通過調(diào)用ArcPy制圖模塊的函數(shù)和類來自動(dòng)完成[3]。

（1）裁剪數(shù)據(jù)源

由于原始Shapefile中包含的是亞洲和俄羅斯的火點(diǎn)數(shù)據(jù)，所以需要從大量的數(shù)據(jù)中篩選出研究區(qū)的火點(diǎn)數(shù)據(jù)進(jìn)行批量處理[4]。首先以哈爾濱的矢量邊界基準(zhǔn)，裁剪掉哈爾濱范圍外的火點(diǎn)數(shù)據(jù)；再利用耕地的矢量邊界為基準(zhǔn)，裁剪掉耕地范圍外的火點(diǎn)數(shù)據(jù)。以下代碼部分以裁剪掉哈爾濱范圍外的火點(diǎn)數(shù)據(jù)為例：

①將哈爾濱邊界的矢量數(shù)據(jù)拷貝到與代碼同一路徑下，并獲取其路徑：

_clipShp = os.getcwd() + "\" + "border" + "\" +"border.shp"

②獲取當(dāng)前地理處理的工作路徑并定義被裁減的數(shù)據(jù)路徑和裁剪后輸出的路徑：

env.workspace = _pathNewDate

inputShp=_pathNewDate+"_"+"24h"+".shp"

outShp=_pathNewDate+"_"+"24h"+"_out.shp"

arcpy.Clip_analysis(inputShp, _clipShp, outShp)

（2）修復(fù)數(shù)據(jù)鏈接并制圖

打開一個(gè)長時(shí)間沒有使用的地圖文檔，在通過ArcMap打開該地圖文檔時(shí)，ArcMap內(nèi)容列表中的圖層就會(huì)有一個(gè)帶著紅色感嘆號(hào)的標(biāo)記，并且數(shù)據(jù)框中不會(huì)顯示任何圖層。

出現(xiàn)這種情況是因?yàn)樵瓉淼臄?shù)據(jù)鏈接已經(jīng)斷開，而造成數(shù)據(jù)鏈接斷開的原因可能有以下幾種：地圖文檔是以絕對(duì)路徑保存，但是數(shù)據(jù)源的路徑已經(jīng)發(fā)生改變。例如，數(shù)據(jù)的源文件被移動(dòng)到其他地方；地圖文檔以相對(duì)路徑保存，但是.mxd文檔被移動(dòng)到其他地方，而數(shù)據(jù)的源文件位置沒有改變；數(shù)據(jù)源的名稱被修改過。

圖2 解壓后的文件在計(jì)算機(jī)中的排列方式

數(shù)據(jù)鏈接斷開的情況經(jīng)常發(fā)生，而且手動(dòng)修復(fù)數(shù)據(jù)鏈接也相當(dāng)繁瑣。一旦確定是哪種原因造成數(shù)據(jù)鏈接的斷開，就可以使用腳本來自動(dòng)修復(fù)數(shù)據(jù)鏈接。腳本可以在不打開地圖文檔的情況下檢測(cè)和修復(fù)斷開的數(shù)據(jù)鏈接。本文中，主要是針對(duì)火點(diǎn)的數(shù)據(jù)鏈接進(jìn)行修復(fù)，由于腳本運(yùn)行時(shí)，每一次火點(diǎn)的數(shù)據(jù)源鏈接都在發(fā)生變化，所以下面Python代碼中，主要是采用替換數(shù)據(jù)源的方式來進(jìn)行修復(fù)工作。

①獲得地圖文檔對(duì)象，并檢測(cè)到斷開的數(shù)據(jù)鏈接的圖層，并進(jìn)行修復(fù)：

_pathMxd = os.getcwd() + "\" + "mxd" + "\" +"baseMap.mxd"

mxd = arcpy.mapping.MapDocument(_pathMxd)

brkLyrList = arcpy.mapping.ListBrokenDataSources(mxd)

fireLyr=brkLyrList[0]

fireLyr.replaceDataSource(_pathNewDate,"SHAPEFILE_WORKSPACE",_date+"_"+"24h"+"_out.shp")

②統(tǒng)計(jì)當(dāng)前研究區(qū)的火點(diǎn)個(gè)數(shù)：

iCount=0

cursor = arcpy.da.SearchCursor(_date+"_"+"24h"+"_out.shp", ['LATITUDE'])

for row in cursor:

iCount=iCount+1

③制作火點(diǎn)專題圖，并輸出到指定文件目錄：

_pathJpg = new_path + "\"+jpg+"\"

title = arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT")[1]

title.text=datetime.datetime.now().strftime("%Y-%m-%d")+"(24小時(shí))"

fireCount=arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT")[0]

fireCount.text=str(iCount)+"個(gè)火點(diǎn)"

_pathNewDateJpg=_pathJpg+_date+"\"

arcpy.mapping.Export ToJPEG(m x d,_pathNewDateJpg+_date+"_"+"24h"+".jpg")

上述代碼運(yùn)行，運(yùn)行結(jié)果產(chǎn)生出截至2018年3月31日24小時(shí)內(nèi)和2018年4月1日的哈爾濱疑似火點(diǎn)分布圖（圖3）。

3.逆地址解析及數(shù)據(jù)存儲(chǔ)

（1）逆地址解析

①由2可知，Shapefile數(shù)據(jù)的屬性數(shù)據(jù)包含經(jīng)度（Latitude），緯度（Longitude），火點(diǎn)的發(fā)生日期（Acq_Date）和時(shí)間（Acq_Time）。在腳本中，使用游標(biāo)訪問經(jīng)過地理處理的矢量數(shù)據(jù)的屬性列表中指定字段的數(shù)據(jù)。ArcPy.da模塊提供了三種類型的游標(biāo)，分別是搜索、插入和更新，在此過程中，程序中主要使用的是搜索游標(biāo)。

cursor = arcpy.da.SearchCursor(_date+"_"+"24h"+"_out.shp", "LATITUDE","LONGITUDE","ACQ_DATE","ACQ_TIME"])

②將上一步驟中搜索到的所有數(shù)據(jù)進(jìn)行for循環(huán)，獲取對(duì)應(yīng)火點(diǎn)的經(jīng)度和緯度：

for row in cursor:

x = row[0]

y = row[1]

ACQ_DATE = row[2]

ACQ_TIME = row[3]

其中，x代表的經(jīng)度，y代表的是緯度。

③調(diào)用騰訊逆地址解析的 API接口，解析每一個(gè)火點(diǎn)的具體地址。該API接口的請(qǐng)求方式是GET請(qǐng)求，其中URL中的主要參數(shù)說明見表1。

圖3 疑似火點(diǎn)分布圖

導(dǎo)入第三方庫requests，它是python實(shí)現(xiàn)的簡(jiǎn)單易用的HTTP庫，requests.get( )用于請(qǐng)求目標(biāo)網(wǎng)站，類型是一個(gè)HTTPresponse類型，返回是一個(gè)具有JSON數(shù)據(jù)格式的字符串。為了更加方便的從字符串中取出地址信息，在腳本中引入JSON模塊，將字符串轉(zhuǎn)換成python中的字典格式的數(shù)據(jù)：

import requests,json

url = "http://apis.map.qq.com/ws/geocoder/v1/?location="+ str(y) + ","+ str(x)+"&key=KGRBZ-CJ2WP-BKYD5-VGPJI-XELGT-5IB23&get_poi=0"

response = requests.get(url).content

results = json.loads(response)

④從字典中取出詳細(xì)的地址信息，由于某些經(jīng)緯度地址無法精確到鄉(xiāng)、鎮(zhèn)、街道或者村、屯，所以在腳本中需要添加判斷語句，防止信息獲取不完整。如果字典的鍵值中包含關(guān)鍵字"town"和"landmark_l2"，則說明該經(jīng)緯度能解析到鄉(xiāng)、鎮(zhèn)、街道，反之，則將Town變量和設(shè)置為占位符號(hào)" - "。村、屯信息的獲取則采用正則表達(dá)式驗(yàn)證字典中的數(shù)據(jù)是否符合要求：

表1 URL主要參數(shù)

P r o v i n c e =results["result"]["address_component"] ["province"]

City = results["result"]["address_component"]["city"]

District = results["result"]["address_component"]["district"]

town = "town”

if (town in results["result"]["address_reference"].keys

()):

Town = results["result"]["address_reference"]["town" ]["title"]

else :

Town = " - "

PATTERN2 = r'([u4e00-u9fa5](?:村|屯))'

landmark_l2 = "landmark_l2"

if (landmark_l2 in results["result"] ["address_reference"].keys()):

Village = results["result"]["address_reference"]["landmark_l2"]["title"].encode("utf-8")

Village = Vil.decode('utf-8')

pattern2 = re.compile(PATTERN2)

m2 = pattern2.search(Vil)

if m2:

V = Village

else:

V = " - "

else:

V = " - "

（2）非結(jié)構(gòu)化數(shù)據(jù)存儲(chǔ)

①在MySQL數(shù)據(jù)庫中新建數(shù)據(jù)庫，在數(shù)據(jù)庫下創(chuàng)建數(shù)據(jù)表表名設(shè)為firedata，數(shù)據(jù)庫結(jié)構(gòu)見表2。

②導(dǎo)入Python第三方庫MySQLdb連接MySQL數(shù)據(jù)庫并獲取相應(yīng)的數(shù)據(jù)表的游標(biāo)，準(zhǔn)備存儲(chǔ)數(shù)據(jù)，并在所有數(shù)據(jù)存儲(chǔ)完成關(guān)閉游標(biāo)和數(shù)據(jù)庫連接[5]，防止腳本運(yùn)行中出現(xiàn)不必要的錯(cuò)誤：

conn = MySQLdb.connect(host="localhost",user="root",passwd="mypassword",charset="gbk")

curs = conn.cursor()

conn.select_db('firedatabase')

tableName="firedata"

value = [x, y, Province, City, District, Town, V,ACQ_DATE，ACQ_TIME]

表2 數(shù)據(jù)庫結(jié)構(gòu)

curs.execut

e("insertinto "+table Name+"(X,Y,P,C,D,T,V,ACQ_DATE,ACQ_TIME) values(%s,%s,%s,%s,%s,%s,%s,%s,%s)" ,value)

conn.commit()

curs.close()

conn.close()

其中，host表示數(shù)據(jù)的主機(jī)名稱或者IP地址；user表示用戶名稱；passwd表示數(shù)據(jù)庫的密碼；charset表示存入數(shù)據(jù)庫的數(shù)據(jù)編碼格式，gbk代表漢字內(nèi)碼擴(kuò)展規(guī)范。

5.結(jié)束語

通過對(duì)腳本運(yùn)行測(cè)試的結(jié)果來分析，為火點(diǎn)解譯工作人員節(jié)省了大量的時(shí)間，而且在數(shù)據(jù)自動(dòng)存儲(chǔ)方面有利于保存歷史數(shù)據(jù)，以便以后的研究人員對(duì)火點(diǎn)產(chǎn)生的原因及時(shí)間段進(jìn)行研究和分析。