shape files convert to rasters

Table of contents

背景
withinD5.py程式設計重點
ToDo
codes
Resource
- links
- notes

背景

目標

本次作業的目標在於將內政部(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個步驟，似過於繁瑣。此處予以整併。

檔案IO及格式說明

shape 檔：（輸入）TOWN_MOI_1090727.shp 來自MOI官網
nc檔：（輸入及輸出）模版template_d5_1x1.nc及最終結果(同名)
csv檔：（輸出）鄉鎮區屬性內容

withinD5.py程式設計重點

shape檔屬性內容之讀取與紀錄

基本上shape檔的骨幹也是以dict型式儲存多邊形座標點，其屬性內容儲存在fields內，如下所示。
以MOI檔案而言，共有7個fields，將其作成DataFrame的columns(col)，
- 每個鄉鎮區紀錄(rec)前7項記住了各區的名稱代碼等。
- 雖然是dict，但shape檔案的fields紀錄仍然是按照順序排列的。

shp='TOWN_MOI_1090727.shp'
shape = shapefile.Reader(shp)
f=shape.fields
rec=shape.records()
col=[i[0] for i in f[1:]]
df=DataFrame({col[ic]:[rec[i][ic] for i in range(len(rec))] for ic in range(7)})
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。

plgs,multi=[],[] #(lon,lat)
for i in range(len(df)):
  tp=shape.shapeRecords()[i].shape.__geo_interface__['type']
  cr=shape.shapeRecords()[i].shape.__geo_interface__['coordinates']
  if len(cr)!=1:
    multi.append(i)
    if type(cr[0][0][0])==tuple:
      crs=[cr[ic][0][:] for ic in range(len(cr))]
    else:
      crs=[cr[ic][:] for ic in range(len(cr))]
    plg=[]
    for cr in crs:
      plg+=cr
    plgs.append(plg)
  else:
    plgs.append(cr[0])

經緯度轉置與邊界範圍

MOI座標是(經度、緯度)，shaple約定是(緯度、經度)，順序相反

Dplg=[] #(lat,lon)
for plg in plgs:
  lon=[i[0] for i in plg[:]]
  lat=[i[1] for i in plg[:]]
  crd=[(i,j) for i,j in zip(lat,lon)]
  Dplg.append(crd)

記錄多邊形邊界經緯度的範圍

    58  mnLat,mnLon,mxLat,mxLon=[np.array([Polygon(plg).bounds[i] for plg in Dplg]) for i in range(4)]

單一多邊形within之判別

對每一個網格點進行判別，共有(414, 252)點，大多為海上點。
由於鄉鎮區總數368，很多是在多邊形範圍外，因此必須進行座標的篩選，
經測試，np.where會比within更快
- 此處做緯度、經度2次np.where判別
- 如果n屬於multi多個多邊形個案，必須每一個多邊形判別，不在此處理。

nplgs=len(df)
DIS=np.zeros(shape=(nrow2,ncol2))
for j in range(nrow2):
  for i in range(ncol2):
    p1=Point((Plat[j,i],Plon[j,i]))
    idx=np.where((Plat[j,i]-mnLat)*(Plat[j,i]-mxLat)<=0)
    if len(idx[0])==0: continue
    idx2=np.where((Plon[j,i]-mnLon[idx[0][:]])*(Plon[j,i]-mxLon[idx[0][:]])<=0)
    if len(idx2[0])==0: continue
    for n in list(idx[0][idx2[0]]):                                              #loop for each polygons
      if n in multi:continue
      poly = Polygon(Dplg[n])
      if p1.within(poly):        #boolean to check whether the p1 coordinates is inside the polygon or not
        DIS[j,i]=float(n)
        break

多個多邊形within之判別

此處以鄉鎮區為主要的迴圈，針對網格點進行篩選，符合邊界範圍條件網格才進行within判別
- 3層迴圈順序、由外而內分別為：鄉鎮區->多邊形->網格點
- 因前述plgs將所有多邊形lump在一起，within判別會出錯，必須重新讀取shape檔、依序判別同一鄉鎮區的多個多邊形。
- 同樣進行層次的確認（78~81行）、與經緯度的轉置(83~85行)

for n in multi:
  idx=np.where((Plat-mnLat[n])*(Plat-mxLat[n])<=0)
  idx2=np.where((Plon[idx[0][:],idx[1][:]]-mnLon[n])*(Plon[idx[0][:],idx[1][:]]-mxLon[n])<=0)
  cr=shape.shapeRecords()[n].shape.__geo_interface__['coordinates']
  if type(cr[0][0][0])==tuple:
    plgs=[cr[ic][0][:] for ic in range(len(cr))]
  else:
    plgs=[cr[ic][:] for ic in range(len(cr))]
  for plg in plgs:
    lon=[i[0] for i in plg[:]]
    lat=[i[1] for i in plg[:]]
    crd=[(ii,jj) for ii,jj in zip(lat,lon)]
    for ij in range(len(idx2[0])):
      j=idx[0][idx2[0][ij]]
      i=idx[1][idx2[0][ij]]
      p1=Point((Plat[j,i],Plon[j,i]))
      poly = Polygon(crd)
      if p1.within(poly):  #boolean to check whether the p1 coordinates is inside the polygon or not
        DIS[j,i]=float(n)

模版製作並輸出檔案

由於nc檔案只能變動unlimited dimension，因此必須先將rec_dmn設成ROW及COL，
第一次要在各維度依序展開(107~110行)
- 以後只要維度相同，可以一次倒入數據，但是bound要設好(112行)
- 只要一個變數及TFLAGS，其餘變數不必留存
- 使用ncrename程式將變數NO成NUM_TOWN

#ncks -O --mk_rec_dmn ROW template_d5_1x1.nc a.nc
#ncks -O --mk_rec_dmn COL b.nc c.nc
#ncks -O -v NO,TFLAG -d TSTEP,0 c.nc template_d5_1x1.nc
#ncrename -v NO,NUM_TOWN $nc
nc = netCDF4.Dataset('template_d5_1x1.nc', 'r+')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
nt,nlay,nrow,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4))
nc.NCOLS=ncol2
nc.NROWS=nrow2
nc.NVARS=1
nc.NSTEPS=1
nc.XCELL=RESm
nc.YCELL=RESm
if nrow!=nrow2 or ncol!=ncol2:
  for j in range(nrow2):
    for i in range(ncol2):
      nc.variables['NUM_TOWN'][0,0,j,i]=DIS[j,i]
else:
  nc.variables['NUM_TOWN'][0,0,:nrow,:ncol]=DIS[:,:]
nc.close()

ToDo

空間模版改變時注意事項
- 中心點位置：(YCENT及XCENT)程式（line 9）與nc檔案都會有，應以後者為主，須調整line 9與line 13讀取的順序
- 此處以模版範圍為主，沒有其他定義
- 網格數：由程式自行計算，也無需另外定義。
- 解析度：@line 20~23, 105~106
shape檔結構、格式的Try and Err
- 同一序號，多個多邊形。不見得其他的shape檔也是此一結構。
- 多個多邊形、不同序列層次套疊。

shape files convert to rasters

背景

目標

方案考量

檔案IO及格式說明

withinD5.py程式設計重點

shape檔屬性內容之讀取與紀錄

多邊形個數之區分

經緯度轉置與邊界範圍

單一多邊形within之判別

多個多邊形within之判別

模版製作並輸出檔案

ToDo

codes

Resource

links

notes