周鳳杰，萬魯河*

(1.哈爾濱師范大學(xué)寒區(qū)地理環(huán)境監(jiān)測(cè)與地理信息服務(wù)黑龍江省重點(diǎn)實(shí)驗(yàn)室；2.哈爾濱師范大學(xué))

0 引言

焚燒秸稈會(huì)在短時(shí)間內(nèi)嚴(yán)重影響到空氣質(zhì)量，污染城市大氣環(huán)境.2017年10月9日，哈爾濱市人民政府印發(fā)《哈爾濱市禁止野外焚燒秸稈改善大氣環(huán)境質(zhì)量實(shí)施方案》(以下簡(jiǎn)稱《方案》)，根據(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).該文通過編寫Python腳本，快速、實(shí)時(shí)的從MODIS衛(wèi)星上獲取火點(diǎn)數(shù)據(jù)，并自動(dòng)實(shí)時(shí)生成火點(diǎn)分布圖，從而更好地解決這一問題.

1 研究區(qū)概況

哈爾濱市處于中國(guó)東北地區(qū)、東北亞中心地帶、黑龍江省南部，地跨東經(jīng)125°42′～130°10′、北緯44°04′～46°40′,地勢(shì)東南高,西北低.根據(jù)2018年哈爾濱市統(tǒng)計(jì)年鑒，2018年年末哈爾濱市占地面積為5.31萬 km2，耕地面積為1.97756萬 km2，約占總面積的37.24%.其中市區(qū)耕地面積為0.540434萬 km2，市轄縣(市)耕地面積為1.437131萬 km2.如圖1所示.

2 數(shù)據(jù)來源

所用的火點(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 h內(nèi)的火點(diǎn)數(shù)據(jù)，48 h的火點(diǎn)數(shù)據(jù)和最近7 d內(nèi)的火點(diǎn)數(shù)據(jù).其中，Google Earth KML，KML 全稱是Keyhole Markup Language KML,是一個(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è)文件：(1)主文件(*.shp)，存儲(chǔ)地理要素的幾何圖形的文件；(2)索引文件(*.shx),存儲(chǔ)圖形要素與屬性信息索引的文件；(3)dBASE表文件(*.dbf)，-存儲(chǔ)要素信息屬性的dBase表文件[1].不同的數(shù)據(jù)格式，但是都用相同的字段標(biāo)識(shí)，文件中共包含13個(gè)字段Latitude，Longitude，Bright_ti4，Scan，Track，Acq_Date，Acq_Time，Satellite，Confidence，Version，Bright_ti5，F(xiàn)RP，DayNight；分別表示火點(diǎn)的經(jīng)度，緯度，亮度溫度I-4，沿掃描像素大小，沿著軌道像素大小，采集日期，采集時(shí)間，衛(wèi)星，置信度，版本(收集和來源，其中“1.0NRT”表示收集1 NRT處理，“1.0” - 集合1標(biāo)準(zhǔn)處理)亮度溫度I-5，火輻射能量，白天或晚上(其中，D表示白天火災(zāi)，N表示夜間火災(zāi)).該文以爬取24 h內(nèi)的Shapefile數(shù)據(jù)為例.

圖1 哈爾濱市耕地分布示意圖

3 Python簡(jiǎn)介

Python是一門解釋型語言，使用專門的解釋器對(duì)源程序進(jìn)行逐行解釋成某個(gè)特定平臺(tái)的機(jī)器碼，在運(yùn)行時(shí)將程序翻譯成機(jī)器語言，解釋型語言相當(dāng)于把編譯型語言中的編譯和解釋過程混合在一起同時(shí)完成.它也常被稱作膠水語言，能夠把用其他語言(如C，C++等)制作的各種模塊，引入到腳本中使用[2].

3.1 豐富的開源庫(kù)

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

3.2 可擴(kuò)展性和可嵌入性

Python可以跨操作平臺(tái)運(yùn)行，即Python程序的核心語言和標(biāo)準(zhǔn)庫(kù)可以在Linux、Windows及其他帶有Python解釋器的平臺(tái)上無差別的運(yùn)行[4].Python的標(biāo)準(zhǔn)實(shí)現(xiàn)是由可移植的ANSIC編寫的.

3.3 代碼簡(jiǎn)單易懂，可讀性強(qiáng)

Python語法簡(jiǎn)潔清晰，突出特點(diǎn)之一就是強(qiáng)制用空白符(white space)作為語句縮進(jìn).相對(duì)于其他高度結(jié)構(gòu)化的編程語言(C++或者Visual Basic)而言，Python更容易被掌握.它的語法簡(jiǎn)單，編程者將有更多的時(shí)間來解決實(shí)際問題，而不需要再學(xué)習(xí)Python語言上花費(fèi)太多的時(shí)間和精力.閱讀一個(gè)良好的Python程序就感覺像是在讀英語一樣，盡管這個(gè)英語的要求非常嚴(yán)格.Python的這種偽代碼本質(zhì)是它最大的優(yōu)點(diǎn)之一[5].

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

4.1 原始數(shù)據(jù)處理

4.1.1 數(shù)據(jù)下載

Python標(biāo)準(zhǔn)庫(kù)中的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)行讀取[6].文中利用urllib2模塊下載Shapefile數(shù)據(jù)格式的壓縮包數(shù)據(jù)，并設(shè)置數(shù)據(jù)存儲(chǔ)的位置.導(dǎo)入urllib2庫(kù)下載壓縮包，并導(dǎo)入Python標(biāo)準(zhǔn)庫(kù)中的os模塊，用于程序中切換文件路徑、自動(dòng)創(chuàng)建或者移除文件夾、更改文件的名字、判斷指定文件夾在指定路徑下是否存在等操作.

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

(1)獲取系統(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)

(2)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) C hrome/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ī)中的排列方式如圖2所示.

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

4.1.2 數(shù)據(jù)解壓以及重命名

導(dǎo)入解壓zip數(shù)據(jù)格式的第三方Python庫(kù)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：

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

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

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

(2)判斷文件路徑是否已經(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"

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

os.chdir(_pathNewDate)

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

newShpDir = os.listdir(_pathNewDate)

(4)利用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ī)中的排列方式如圖3所示.

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

4.2 地理數(shù)據(jù)處理

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

4.2.1 裁剪數(shù)據(jù)源

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

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

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

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

env.workspace = _pathNewDate

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

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

arcpy.Clip_analysis(inputShp, _clipShp, outShp)

4.2.2 修復(fù)數(shù)據(jù)鏈接并制圖

打開一個(gè)長(zhǎng)時(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ù)源的名稱被修改過.

數(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ù)工作.

(1)獲得地圖文檔對(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")

(2)統(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

(3)制作火點(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.ExportToJPEG(mxd, _pathNewDateJpg+

_date+"_"+"24h"+".jpg")

上述代碼運(yùn)行，運(yùn)行結(jié)果產(chǎn)生的專題圖如圖4所示，表示的是截止2018年3月31日24 h內(nèi)(左)和2018年4月1日(右)的哈爾濱市疑似火點(diǎn)分布圖.

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

4.3.1 逆地址解析

(1)由4.2可知，Shapefile數(shù)據(jù)的屬性數(shù)據(jù)包含經(jīng)度(Latitude)，緯度(Longitude)，火點(diǎn)的發(fā)生日期(Acq_Date)和時(shí)間(Acq_Time).在腳本中，使用游標(biāo)訪問4.2中經(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"])

(2)將上一步驟中搜索到的所有數(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]

圖4 哈爾濱市疑似火點(diǎn)分布示意圖

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

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

表1 URL的主要參數(shù)說明

導(dǎo)入第三方庫(kù)requests,它是Python實(shí)現(xiàn)的簡(jiǎn)單易用的HTTP庫(kù)，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)

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

Province = 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 = " - "

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

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

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

(2)導(dǎo)入Python第三方庫(kù)MySQLdb連接MySQL數(shù)據(jù)庫(kù)并獲取相應(yīng)的數(shù)據(jù)表的游標(biāo)，準(zhǔn)備存儲(chǔ)數(shù)據(jù)，并在所有數(shù)據(jù)存儲(chǔ)完成(4.3.1中for循環(huán)結(jié)束)關(guān)閉游標(biāo)和數(shù)據(jù)庫(kù)連接[10]，防止腳本運(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]

curs.execute("insert into "+tableName+"(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ù)庫(kù)的密碼；charset表示存入數(shù)據(jù)庫(kù)的數(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)行研究和分析.