Table of contents 本次作業的目標在於將內政部(MOI 鄉鎮區界的shape檔，讀成d5範圍1公里 解析度的raster檔(ioapi nc格式)備用。 shape file的讀取方式有很多種(詳參GIS Stack Exchange )，大多藉由gdal程式庫 。然而因gdal套用其他軟體相依性太高，太過複雜，裝設不易，且有記憶體容量之限制。 一般寫出方法以tiff為主，目前尚沒有nc的現成套件。 此處選擇直讀式的PyShp ，其裝置、使用範例可以參考PyPi官網 。PyShp 的弱點在沒有現成的寫出方法可以套用，必須另外自行撰寫。然因作業所要求格式對GIS領域而言過於冷僻，因此勢必無法規避自行撰寫。 vector to raster 的核心是網格點是否在多邊形之內的判定(within)，可以調用shapely 的程式庫，詳官網說明及及python解析KML(GML)檔 範例。 kml也可以轉成raster，一般shp向量檔也有kml形式儲存備用的。然而kml的轉檔需執行rd_kml 及withinD1 等檔案讀取與寫出2個步驟，似過於繁瑣。此處予以整併。 shape 檔：（輸入）TOWN_MOI_1090727.shp 來自MOI官網 nc檔：（輸入及輸出）模版template_d5_1x1.nc及最終結果(同名)csv檔：（輸出）鄉鎮區屬性內容 基本上shape檔的骨幹也是以dict型式儲存多邊形座標點，其屬性內容儲存在fields內，如下所示。 以MOI 檔案而言，共有7個fields，將其作成DataFrame的columns(col)，每個鄉鎮區紀錄(rec)前7項記住了各區的名稱代碼等。 雖然是dict，但shape檔案的fields紀錄仍然是按照順序排列的。 28 shp = 'TOWN_MOI_1090727.shp'
29 shape = shapefile . Reader ( shp )
30 f = shape . fields
31 rec = shape . records ()
32 col = [ i [ 0 ] for i in f [ 1 :]]
33 df = DataFrame ({ col [ ic ]:[ rec [ i ][ ic ] for i in range ( len ( rec ))] for ic in range ( 7 )})
34 df . to_csv ( 'record.csv' )
個數 之區分 由於鄉鎮區範圍可能是接壤的單一多邊形 、或者不連續、多個多邊形 ，單一多邊形 處理較為單純且佔多數，多個多邊形 較為複雜但個數較少。在迴圈的設計上分別處理。先記錄前者(調用shape.__geo_interface__['coordinates'])及multi標籤。 plgs順序必須保持跟df相同，multi則不必。 在multi之24個個案中，MOI 格式略有不同，原因不明。經嘗試錯誤發現，有的是2層套疊list，有的是單層，必須探究哪一層才是經緯度tuple，以此來判別。 先將多個多邊形lump成一個序列(第46~48行)、保持正確的格式、能夠計算Polygon.bounds就好，後面會再依照多邊形個數依序讀取、判別within。 36 plgs , multi = [],[] #(lon,lat)
37 for i in range ( len ( df )):
38 tp = shape . shapeRecords ()[ i ]. shape . __geo_interface__ [ 'type' ]
39 cr = shape . shapeRecords ()[ i ]. shape . __geo_interface__ [ 'coordinates' ]
40 if len ( cr ) != 1 :
41 multi . append ( i )
42 if type ( cr [ 0 ][ 0 ][ 0 ]) == tuple :
43 crs = [ cr [ ic ][ 0 ][:] for ic in range ( len ( cr ))]
44 else :
45 crs = [ cr [ ic ][:] for ic in range ( len ( cr ))]
46 plg = []
47 for cr in crs :
48 plg += cr
49 plgs . append ( plg )
50 else :
51 plgs . append ( cr [ 0 ])
MOI 座標是(經度、緯度)，shaple約定是(緯度、經度)，順序相反 52 Dplg = [] #(lat,lon)
53 for plg in plgs :
54 lon = [ i [ 0 ] for i in plg [:]]
55 lat = [ i [ 1 ] for i in plg [:]]
56 crd = [( i , j ) for i , j in zip ( lat , lon )]
57 Dplg . append ( crd )
58 mnLat , mnLon , mxLat , mxLon = [ np . array ([ Polygon ( plg ). bounds [ i ] for plg in Dplg ]) for i in range ( 4 )]
對每一個網格點進行判別，共有(414, 252)點，大多為海上點。 由於鄉鎮區總數368，很多是在多邊形範圍外，因此必須進行座標的篩選， 經測試，np.where會比within更快此處做緯度、經度2次np.where判別 如果n屬於multi多個多邊形個案，必須每一個多邊形判別，不在此處理。 59 nplgs = len ( df )
60 DIS = np . zeros ( shape = ( nrow2 , ncol2 ))
61 for j in range ( nrow2 ):
62 for i in range ( ncol2 ):
63 p1 = Point (( Plat [ j , i ], Plon [ j , i ]))
64 idx = np . where (( Plat [ j , i ] - mnLat ) * ( Plat [ j , i ] - mxLat ) <= 0 )
65 if len ( idx [ 0 ]) == 0 : continue
66 idx2 = np . where (( Plon [ j , i ] - mnLon [ idx [ 0 ][:]]) * ( Plon [ j , i ] - mxLon [ idx [ 0 ][:]]) <= 0 )
67 if len ( idx2 [ 0 ]) == 0 : continue
68 for n in list ( idx [ 0 ][ idx2 [ 0 ]]): #loop for each polygons
69 if n in multi : continue
70 poly = Polygon ( Dplg [ n ])
71 if p1 . within ( poly ): #boolean to check whether the p1 coordinates is inside the polygon or not
72 DIS [ j , i ] = float ( n )
73 break
此處以鄉鎮區為主要的迴圈，針對網格點進行篩選，符合邊界範圍條件網格才進行within判別3層迴圈順序、由外而內分別為：鄉鎮區->多邊形->網格點 因前述plgs將所有多邊形lump在一起，within判別會出錯，必須重新讀取shape檔、依序判別同一鄉鎮區的多個多邊形。 同樣進行層次的確認（78~81行）、與經緯度的轉置(83~85行) 74 for n in multi :
75 idx = np . where (( Plat - mnLat [ n ]) * ( Plat - mxLat [ n ]) <= 0 )
76 idx2 = np . where (( Plon [ idx [ 0 ][:], idx [ 1 ][:]] - mnLon [ n ]) * ( Plon [ idx [ 0 ][:], idx [ 1 ][:]] - mxLon [ n ]) <= 0 )
77 cr = shape . shapeRecords ()[ n ]. shape . __geo_interface__ [ 'coordinates' ]
78 if type ( cr [ 0 ][ 0 ][ 0 ]) == tuple :
79 plgs = [ cr [ ic ][ 0 ][:] for ic in range ( len ( cr ))]
80 else :
81 plgs = [ cr [ ic ][:] for ic in range ( len ( cr ))]
82 for plg in plgs :
83 lon = [ i [ 0 ] for i in plg [:]]
84 lat = [ i [ 1 ] for i in plg [:]]
85 crd = [( ii , jj ) for ii , jj in zip ( lat , lon )]
86 for ij in range ( len ( idx2 [ 0 ])):
87 j = idx [ 0 ][ idx2 [ 0 ][ ij ]]
88 i = idx [ 1 ][ idx2 [ 0 ][ ij ]]
89 p1 = Point (( Plat [ j , i ], Plon [ j , i ]))
90 poly = Polygon ( crd )
91 if p1 . within ( poly ): #boolean to check whether the p1 coordinates is inside the polygon or not
92 DIS [ j , i ] = float ( n )
由於nc檔案只能變動unlimited dimension，因此必須先將rec_dmn設成ROW及COL， 第一次要在各維度依序展開(107~110行)以後只要維度相同，可以一次倒入數據，但是bound要設好(112行) 只要一個變數及TFLAGS，其餘變數不必留存 使用ncrename程式將變數NO成NUM_TOWN 94 #ncks -O --mk_rec_dmn ROW template_d5_1x1.nc a.nc
95 #ncks -O --mk_rec_dmn COL b.nc c.nc
96 #ncks -O -v NO,TFLAG -d TSTEP,0 c.nc template_d5_1x1.nc
97 #ncrename -v NO,NUM_TOWN $nc
98 nc = netCDF4 . Dataset ( 'template_d5_1x1.nc' , 'r+' )
99 V = [ list ( filter ( lambda x : nc . variables [ x ]. ndim == j , [ i for i in nc . variables ])) for j in [ 1 , 2 , 3 , 4 ]]
100 nt , nlay , nrow , ncol = ( nc . variables [ V [ 3 ][ 0 ]]. shape [ i ] for i in range ( 4 ))
101 nc . NCOLS = ncol2
102 nc . NROWS = nrow2
103 nc . NVARS = 1
104 nc . NSTEPS = 1
105 nc . XCELL = RESm
106 nc . YCELL = RESm
107 if nrow != nrow2 or ncol != ncol2 :
108 for j in range ( nrow2 ):
109 for i in range ( ncol2 ):
110 nc . variables [ 'NUM_TOWN' ][ 0 , 0 , j , i ] = DIS [ j , i ]
111 else :
112 nc . variables [ 'NUM_TOWN' ][ 0 , 0 ,: nrow ,: ncol ] = DIS [:,:]
113 nc . close ()
空間模版改變時注意事項中心點位置：(YCENT及XCENT)程式（line 9）與nc檔案都會有，應以後者為主，須調整line 9與line 13讀取的順序 此處以模版範圍為主，沒有其他定義 網格數：由程式自行計算，也無需另外定義。 解析度：@line 20~23, 105~106 shape檔結構、格式的Try and Err同一序號，多個多邊形。不見得其他的shape檔也是此一結構。 多個多邊形、不同序列層次套疊。 disscussion, How to read a shapefile in Python? , Stack Exchange , Sep 28 2018. jlawhead, PyShp 2.1.3 , pypi , Jan 14, 2021 Henrikki Tenkanen, Point in Polygon & Intersect , github.io , Jan 17, 2018. Sean Gillies, Shapely 1.8.0 pypi , Oct 26, 2021 Sean Gillies, The Shapely User Manual , readthedocs , Dec 10, 2021