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

周鳳杰，萬魯河*

(1.哈爾濱師范大學寒區(qū)地理環(huán)境監(jiān)測與地理信息服務黑龍江省重點實驗室；2.哈爾濱師范大學)

0 引言

焚燒秸稈會在短時間內嚴重影響到空氣質量，污染城市大氣環(huán)境.2017年10月9日，哈爾濱市人民政府印發(fā)《哈爾濱市禁止野外焚燒秸稈改善大氣環(huán)境質量實施方案》(以下簡稱《方案》)，根據(jù)《方案》，哈爾濱市將建立市級領導包區(qū)縣(市)、區(qū)縣(市)干部包鄉(xiāng)(鎮(zhèn))、鄉(xiāng)(鎮(zhèn))干部包村、鄉(xiāng)(鎮(zhèn))村(屯)干部包農戶和地塊的逐級包保責任制，形成以鄉(xiāng)(鎮(zhèn))、村(屯)、農戶、地塊為單元的網格化管理體系.快速獲取火點數(shù)據(jù)并對秸稈焚燒發(fā)生的地點進行實時核查，成為落實逐級包保責任制關鍵的一環(huán).該文通過編寫Python腳本，快速、實時的從MODIS衛(wèi)星上獲取火點數(shù)據(jù)，并自動實時生成火點分布圖，從而更好地解決這一問題.

1 研究區(qū)概況

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

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

所用的火點經緯度數(shù)據(jù)來自網站https://firms.modaps.eosdis.nasa.gov/，該網站一共提供了3種不同的數(shù)據(jù)格式，包括Google Earth KML,CSV表格，Shapefile；并且每種數(shù)據(jù)格式都提供了3種不同時間段的數(shù)據(jù)，包括24 h內的火點數(shù)據(jù)，48 h的火點數(shù)據(jù)和最近7 d內的火點數(shù)據(jù).其中，Google Earth KML，KML 全稱是Keyhole Markup Language KML,是一個基于XML語法和文件格式的文件，用來描述和保存地理信息如點、線、圖片、折線并在Google Earth 客戶端顯示；CSV(Comma-Separated Values)，逗號分隔值，是一種通用的、相對簡單的文件格式，其文件以純文本形式存儲表格數(shù)據(jù)(數(shù)字和文本)；Shapefile文件是描述控件數(shù)據(jù)的幾何和屬性特征的非拓撲實體矢量數(shù)據(jù)結構的一種格式，一個Shapefile文件至少包括三個文件：(1)主文件(*.shp)，存儲地理要素的幾何圖形的文件；(2)索引文件(*.shx),存儲圖形要素與屬性信息索引的文件；(3)dBASE表文件(*.dbf)，-存儲要素信息屬性的dBase表文件[1].不同的數(shù)據(jù)格式，但是都用相同的字段標識，文件中共包含13個字段Latitude，Longitude，Bright_ti4，Scan，Track，Acq_Date，Acq_Time，Satellite，Confidence，Version，Bright_ti5，F(xiàn)RP，DayNight；分別表示火點的經度，緯度，亮度溫度I-4，沿掃描像素大小，沿著軌道像素大小，采集日期，采集時間，衛(wèi)星，置信度，版本(收集和來源，其中“1.0NRT”表示收集1 NRT處理，“1.0” - 集合1標準處理)亮度溫度I-5，火輻射能量，白天或晚上(其中，D表示白天火災，N表示夜間火災).該文以爬取24 h內的Shapefile數(shù)據(jù)為例.

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

3 Python簡介

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

3.1 豐富的開源庫

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

3.2 可擴展性和可嵌入性

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

3.3 代碼簡單易懂，可讀性強

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

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

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

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

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

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

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

(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變量是一個包含各種瀏覽器信息的字典，目的是用來偽裝瀏覽器爬取數(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)

上述代碼在程序中下載的壓縮包在計算機中的排列方式如圖2所示.

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

4.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：

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

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

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

(2)判斷文件路徑是否已經存在，若不存在，則創(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)切換當前路徑到shp文件存儲的路徑下，并將耕地的空間參考信息prj文件拷貝到路徑下：

os.chdir(_pathNewDate)

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

newShpDir = os.listdir(_pathNewDate)

(4)利用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)

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

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

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

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

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

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

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

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

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

env.workspace = _pathNewDate

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

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

arcpy.Clip_analysis(inputShp, _clipShp, outShp)

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

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

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

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

(1)獲得地圖文檔對象，并檢測到斷開的數(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")

(2)統(tǒng)計當前研究區(qū)的火點個數(shù)：

iCount=0

cursor = arcpy.da.SearchCursor(_date+"_"+

"24h"+"_out.shp", [′LATITUDE′])

for row in cursor:

iCount=iCount+1

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

_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.ExportToJPEG(mxd, _pathNewDateJpg+

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

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

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

4.3.1 逆地址解析

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

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

"ACQ_DATE","ACQ_TIME"])

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

for row in cursor:

x = row[0]

y = row[1]

ACQ_DATE = row[2]

ACQ_TIME = row[3]

圖4 哈爾濱市疑似火點分布示意圖

其中，x代表的經度，y代表的是緯度.

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

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

導入第三方庫requests,它是Python實現(xiàn)的簡單易用的HTTP庫，requests.get()用于請求目標網站，類型是一個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)

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

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

表2 數(shù)據(jù)庫結構說明表

(2)導入Python第三方庫MySQLdb連接MySQL數(shù)據(jù)庫并獲取相應的數(shù)據(jù)表的游標，準備存儲數(shù)據(jù)，并在所有數(shù)據(jù)存儲完成(4.3.1中for循環(huán)結束)關閉游標和數(shù)據(jù)庫連接[10]，防止腳本運行中出現(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]

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

5 結束語

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