基于Python的MODIS衛(wèi)星秸稈焚燒監(jiān)測
——以哈爾濱為例

付嬈

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

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

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

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

所用的火點經(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ù)據(jù)，包括24小時內(nèi)的火點數(shù)據(jù)，48小時的火點數(shù)據(jù)和最近7天內(nèi)的火點數(shù)據(jù)。

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

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

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

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

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

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

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

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

導入Python的標準庫time模塊獲取當前日期與時間，并轉換日期的格式。

①獲取系統(tǒng)時間，并將日期設為全局變量，以日期命名所有下載的文件，以便以后查詢相應時間段的數(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變量是一個包含各種瀏覽器信息的字典，目的是用來偽裝瀏覽器爬取數(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)

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

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

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

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

_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"

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

os.chdir(_pathNewDate)

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

圖1 壓縮包在計算機中的排列方式

newShpDir = os.listdir(_pathNewDate)

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

for temp in newShpDir:

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

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

os.rename(temp,new_name)

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

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

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

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

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

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

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

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

env.workspace = _pathNewDate

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

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

arcpy.Clip_analysis(inputShp, _clipShp, outShp)

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

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

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

圖2 解壓后的文件在計算機中的排列方式

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

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

_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)計當前研究區(qū)的火點個數(shù)：

iCount=0

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

for row in cursor:

iCount=iCount+1

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

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

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

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

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

fireCount.text=str(iCount)+"個火點"

_pathNewDateJpg=_pathJpg+_date+"\"

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

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

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

（1）逆地址解析

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

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

②將上一步驟中搜索到的所有數(shù)據(jù)進行for循環(huá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接口，解析每一個火點的具體地址。該API接口的請求方式是GET請求，其中URL中的主要參數(shù)說明見表1。

圖3 疑似火點分布圖

導入第三方庫requests，它是python實現(xiàn)的簡單易用的HTTP庫，requests.get( )用于請求目標網(wǎng)站，類型是一個HTTPresponse類型，返回是一個具有JSON數(shù)據(jù)格式的字符串。為了更加方便的從字符串中取出地址信息，在腳本中引入JSON模塊，將字符串轉換成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)

④從字典中取出詳細的地址信息，由于某些經(jīng)緯度地址無法精確到鄉(xiāng)、鎮(zhèn)、街道或者村、屯，所以在腳本中需要添加判斷語句，防止信息獲取不完整。如果字典的鍵值中包含關鍵字"town"和"landmark_l2"，則說明該經(jīng)緯度能解析到鄉(xiāng)、鎮(zhèn)、街道，反之，則將Town變量和設置為占位符號" - "。村、屯信息的獲取則采用正則表達式驗證字典中的數(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）非結構化數(shù)據(jù)存儲

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

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

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ù)庫結構

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ù)的主機名稱或者IP地址；user表示用戶名稱；passwd表示數(shù)據(jù)庫的密碼；charset表示存入數(shù)據(jù)庫的數(shù)據(jù)編碼格式，gbk代表漢字內(nèi)碼擴展規(guī)范。

5.結束語

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