shape files convert to rasters

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_kmlwithinD1等檔案讀取與寫出2個步驟,似過於繁瑣。此處予以整併。

檔案IO及格式說明

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

withinD5.py程式設計重點

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

  • 基本上shape檔的骨幹也是以dict型式儲存多邊形座標點,其屬性內容儲存在fields內,如下所示。
  • MOI檔案而言,共有7個fields,將其作成DataFramecolumns(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)]

單一多邊形within之判別

  • 對每一個網格點進行判別,共有(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之判別

  • 此處以鄉鎮區為主要的迴圈,針對網格點進行篩選,符合邊界範圍條件網格才進行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設成ROWCOL
  • 第一次要在各維度依序展開(107~110行)
    • 以後只要維度相同,可以一次倒入數據,但是bound要設好(112行)
    • 只要一個變數及TFLAGS,其餘變數不必留存
    • 使用ncrename程式將變數NONUM_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()

ToDo

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

codes

Resource

  • 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

notes