只在此山中，雲深不知處

聽首歌



© 2018 by Shawn Huang
Last Updated: 2018.5.27

Transtech

QGIS

免費的GIS軟體，首先至官網下載。有多個版本可以選擇，先選擇Standalone installers (MSI) from OSGeo4W packages (recommended for new users)這個，聽說安裝比較簡單，目前版本Version 3.20。以後有機會再試試看Network Installer。下載完成後雙擊檔案安裝，原則上無腦next到底應該即可。有一點要注意的是QGIS已經不支援Win7，所以例如python scripting便不支援，所以作業系統應勿使用Win7，此處為使用Win10。
安裝完成後，找到執行檔開啟。在上方menu bar找到View>Penals & Toolbars來設定快捷鍵的版面設計，以方便工作為準則。如下:

Basics


首先熟悉一下作業環境，先找到下圖位置打開世界地圖，雙擊或是拖曳皆可。
使用View>內的指令移動縮放地圖。可以使用滑鼠滾輪縮放地圖，游標在地圖上移動，下方狀態列的Coordinate會顯示經緯度。
這個地圖原則上就是一張圖，基本上我們可以使用QGIS開啟圖形，開啟後仍然為一個layer，可以在螢幕左下方的Layers處看出。不過隨意開的圖形因為沒有經緯度設計，似乎其左上角座標預設為0,0。

Digital Map


我們需要電子地圖才能操作地圖，光是圖片無法滿足需求。因此我們首先到Natural Earth這個網站下載地圖，選擇SHP或是GeoPackage應該都可以(按滑鼠右鍵下載)。

之後建立一個資料夾解壓縮。為了方便使用，我們可以在左方的Browser處找到Favorites，按滑鼠右鍵選擇Add a Directory，然後選擇地圖位置的資料夾，這個方式可以讓我們較為快速選擇我們需要的檔案。

接下來試著開啟下載的地圖，找到以下檔案開啟:

可以看到世界地圖的國家輪廓，此一圖層僅記錄國家資訊。選擇Layer > Open Attribute Table (F6)或是在圖層上點擊滑鼠右鍵 > Open Attribute Table，可以開啟圖層屬性表格。選擇其中某一或多個資料(使用Shift & Ctrl協助)，回到地圖，使用Zoom to Selection，可以看到被選擇的國家顏色改變，表示被選取。也可以使用快捷列:


Styling


打開Layer > Layer Properties...或是在圖層點擊滑鼠右鍵選擇Properties，然後選擇Symbology。上方顯示Single Symbol，直接改變Color及Opacity等相關屬性。接著將Single Symbol改變為Categorized，Value選擇Continent，然後在下方點擊Classify，出現各不同大陸的顏色列表，可以個別改變顏色，點選套用來改變Style。

再次開啟Properties，選擇Labels，Value選擇NAME，可作修飾後選確定，此時可看到每個國家的名稱。

選擇Layer > Filter...，開啟後輸入如下，點擊確定。此時僅顯示亞洲。

在Layer上點選右鍵 > Export > Save Features As...，出現Save Vector Layer As...視窗，修改如下:

確定後會出現新的圖層。回到原有圖層，開啟Filter...，選擇Clear然後確定，去除原有filter條件，會顯現出世界地圖。將圖層前之勾選去除，此時僅顯現asia圖層。開啟Open Attribute Table，可以看出此圖層的Continent全部都是Asia。

開啟asia圖層之properties > Symbology，與之前之步驟同，根據fid將每個國家著色。此時每個亞洲國家顯示顏色各不相同。在Map_Taiwan資料夾選擇高雄市 > info > EROAD.shp，雙擊或拖曳開啟。

可以在Layers看到開啟EROAD圖層，此時在地圖上的高雄位置應該要有此圖層顯示，不過卻沒看到，不過選擇EROAD圖層然後點Zoom to Layer，應會顯示高雄市街道地圖。開啟Layer > Set CRS of Layer(s)或是右鍵點擊圖層，選Layer CRS > Set Layer CRS...，開啟以下視窗。

Filter輸入TW，選擇合適之經緯度系統，點擊確定，即可在正確位置找到圖層。

讀取檔案資料


首先需有跟地圖相關的資料檔案，因為之後要顯示在地圖，所以最重要的要有x,y位置資料(也就是Longitude & Latitude)，例如:

將其存為tab分隔之文字檔。接著開啟Data Source Manager如下:

Encoding選big5的原因是因為裡面打了中文，注意要對應經緯度。確定後應即可看到對應的都市點。此時可以加上style & label。

Join


設計如下之資料檔案:

將其開啟。選擇cities圖層，右鍵選擇Properties，然後選擇Joins。

點擊下方的+號，出現視窗，選擇要連結的檔案以及對應的欄位，確定後可得到joined view。開啟Attribute Table查看內容，並將人口資料作為Label顯示。
點選cities右鍵 > Export > Save Features As...，將合併後的檔案存為shp file，之後可以直接使用此shp file。

print layout


開啟Porject > New Print Layout...，建立名稱後，進入如下視窗:

選擇Add Item > Add Map加入地圖，然後逐步加入所需的地圖元件。最後可以選擇Layout > Export as Image...，將layout存為圖片。

Selection


Selection可以使用以下按鈕:

或是在Edit > Select > Select Feature(s)，直接框取即可選擇。到Project > Properties... > General > General Settings，可以改變selection的顏色。
在Layer上點選右鍵，選擇Export > Save Selected Features As...

如此可以簡單地取得某一部分的地圖資料。

Merge maps


若欲將兩個file合併為一個，首先開啟此兩個檔案。到Vector > Data Management Tools > Merge vector Layers...開啟如下視窗:

點擊在Input layers後的...，選擇要合併的layer。點擊Merged後的...，決定合併後的檔案類型及名稱。點擊Run執行。執行後便可得到合併在一起的檔案。開啟Attribute Table，觀察若是中文字顯示有問題，點選layer右鍵，開啟properties...，找到Source > Settings，改變encoding(常用的可能為big5或是UTF-8)，應就可正確顯示。

shortest path


計算路網中的最短路徑，首先開啟一個線的圖層。開啟Processing > Toolbox，搜尋Network Analysis，找到Shortest path (point to point)，選擇開啟。


設定network layer， 點擊start point後的...，選擇起點。點擊end point後的...，選擇終點。點擊Run計算，會出現一個暫時layer來顯示路徑。可以修改其style。

PyQGIS

QGIS支援使用Python設計，首先開啟Plugins > Python Console，出現如下之python console:

在>>>後輸入指令，例如 
print("Hello, World")
或是其他簡單的python語法，可正確運作即可。原則上我們常不會只需要一行指令，所以需要編輯器來幫我們寫成一個檔案，QGIS中提供方便的編輯工具，如下:

開啟圖層

因為在QGIS中，處理的原則上是layer的資料，所以首先我們使用程式開啟layer。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
iface.addVectorLayer(uri, "countries", "ogr")
uri儲存的是layer路徑，原則上我們可以在Layer Properties > Information > Information from provider中的Source得到。而iface大概是最重要的物件，用來與QGIS互動。首先使用其中的方法 
addVectorLayer(self, vectorLayerPath: str, baseName: str, providerKey: str) → QgsVectorLayer
參數vectorLayerPath便是uri，接下來的參數baseName是圖層名稱，自選。最後的參數providerKey，對大部分的vector data，通常的值是OGR(也可能是“postgres”, “delimitedtext”, “gpx”, “spatialite”, and “WFS”，可詳見官網說明)。
點擊執行(Run Script)後，應可見圖層被開啟(因為此處執行並不會自動存檔，請記得儲存)。

根據addVectorLayer()的定義，呼叫此方法會傳回QgsVectorLayer物件，也就是圖層，所以我們可以建立以下程式: 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
vlayer = iface.addVectorLayer(uri, "countries", "ogr")
iface.showAttributeTable(vlayer)
方法showAttributeTable()用來開啟Attribute Table。

有的圖檔沒有source，可使用另一個開啟layer的方法，如下: 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
vlayer = iface.addVectorLayer(uri, "countries", "ogr")
ulayer = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
if not ulayer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(ulayer)
iface.showAttributeTable(vlayer)
使用QgsVectorLayer()，變數提供路徑。之後使用QgsProject.instance().addMapLayer(ulayer)語法添加圖層。

我們可以直接開啟圖層，不一定要用程式開啟。此時若要連結到圖層，則使用如下語法: 
ulayer = iface.activeLayer()

印出所有開啟的圖層，使用QgsProject.instance().mapLayers()來取得所有layer，為一個dict。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
vlayer = iface.addVectorLayer(uri, "countries", "ogr")
ulayer = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
if not ulayer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(ulayer)

## 取得所有圖層
layers = QgsProject.instance().mapLayers()
print("----------layers------------")
print(layers)

## 取得layer names
layernames = [lay.name() for lay in QgsProject.instance().mapLayers().values()]
print("----------layernames------------")

## 取得圖層內資料
print(f"{layernames}")
# dictionary with key = layer name and value = layer object
layers_list = {}
for lay in QgsProject.instance().mapLayers().values():
  layers_list[lay.name()] = lay
print("-----------layers_list-----------")
print(layers_list)

# 直接根據名稱取得圖層物件
print("-----------mapLayersByName-----------")
print(QgsProject.instance().mapLayersByName("KH")[0])

所有圖層也是一個樹狀結構，其中的節點有兩種型態:group nodes (QgsLayerTreeGroup) and layer nodes (QgsLayerTreeLayer)。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
vlayer = iface.addVectorLayer(uri, "countries", "ogr")
ulayer = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
if not ulayer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(ulayer)


## 取得root
root = QgsProject.instance().layerTreeRoot()
print("----------root------------")
print(root)
print("----------root.children()------------")
print(root.children())
print(root.children()[0]) # retrieve其中一個child
## 取得圖層id
ids = root.findLayerIds()
print("----------ids------------")
print(ids)
print(ids[0]) # retrieve其中一個id

## 列出Table of content(TOC)中所有checked layers
print("----------checkedLayers------------")
print(root.checkedLayers())

###############################################
## 使用root加入新圖層 
# create a temporary layer
layerPath = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean"
ocean = QgsVectorLayer(layerPath, "Ocean", "ogr")

## 任選一種方式:
# 1. 加到最下方
root.addLayer(ocean)

# 2. 加到特定位置
#root.insertLayer(0, ocean)

# 3. 當然也可以使用如下方式:如前所述
#QgsProject.instance().addMapLayer(ocean)
### -------------------------------------------------- ###

## vector layer與tree layer的轉換
node_layer = root.findLayer(vlayer.id())
print("Layer node:", node_layer) ## << tree node layer
print("Map layer:", node_layer.layer()) ## << map vector layer

## 使用ids來取得所有layer
#ids = root.findLayerIds()
#for i in ids:
#    node_layer = root.findLayer(i)
#    print("Layer node:", node_layer) ## << tree node layer
#    print("Map layer:", node_layer.layer()) ## << map vector layer


## 增加group並將某一layer移至此group
group1 = root.addGroup('One Group') ## 加到最後
#root.insertChildNode(0, group1) ## 加到位置0
theulayer = root.findLayer(ulayer.id()) ## 找到ulayer
cloneu = theulayer.clone() ## clone ulayer
parent = theulayer.parent() # 找到ulayer的parent
group1.insertChildNode(0, cloneu) # 將cloneu加到group1
parent.removeChildNode(theulayer) # 刪除theulayer節點

Attribute Table

此刻的vlayer與ulayer分別代表兩個圖層，我們可以得到其中Attribute Table的資料: 
#uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
#vlayer = iface.addVectorLayer(uri, "countries", "ogr")
#ulayer = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
#if not ulayer.isValid():
#    print("Layer failed to load!")
#else:
#    QgsProject.instance().addMapLayer(ulayer)
iface.showAttributeTable(vlayer)
for field in vlayer.fields():
    print(field.name())
已經完成的步驟可以comment掉，便不會重複執行。valyer.fields()傳回Attribute Table的fields，根據每個field呼叫name()可以得到field的名稱。
類似的觀念，我們可以使用以下語法取得某一欄位的所有資料: 
for feature in vlayer.getFeatures():
    print(feature['fid'])

接下來首先我們先嘗試右鍵點擊layer > filter...，在Query Builder內輸入"ADMIN" like "A%"，此時我們可以看到僅顯示A開頭的國家。選擇Clear清除。接下來使用PyQGIS方法: 
vlayer.setSubsetString("ADMIN like 'A%'")
for feature in vlayer.getFeatures():
    print(feature['ADMIN'])
此時可以看到僅顯示A開頭的國家。回復顯示所有的方式，僅需改為setSubsetString("")即可。

練習: 
for feature in vlayer.getFeatures():
    print(f"{feature['ADMIN']}({feature['NAME_ZH']}) has population = {feature['POP_EST']}")

先練習一下 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

## 屬性欄fields
for field in countries.fields():
    print(field.name(), field.typeName())

print(countries.displayField()) ## 等同於layer properties > Display (點擊Show Map Tips，當滑鼠移至某區塊，將出現該資訊)

取得feature的type。
enum GeometryType { PointGeometry , LineGeometry , PolygonGeometry , UnknownGeometry , NullGeometry }
Note: print x可能導致程式當機，可能是資料量太大 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

########################################################

features = ocean.getFeatures() ## 所有的feature(record)
for feature in features:
#    print("ID:", feature.id(), "Name:", feature['NAME'])
    # 取得每個feature(record)的幾何資料
    geom = feature.geometry()
#    print(geom) ## <QgsGeometry: Polygon(...)>
#    print(geom.type()) ## 2
    geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
#    print(geomSingleType) # True, False
    if geom.type() == QgsWkbTypes.PointGeometry:
            # the geometry type can be of single or multi type
            if geomSingleType: 
                x = geom.asPoint()## shape’s final geometry
                print("Point: ", x)
            else:
                x = geom.asMultiPoint() ## shape’s final geometry
                print("MultiPoint: ", x)
    elif geom.type() == QgsWkbTypes.LineGeometry:
        if geomSingleType:
            x = geom.asPolyline()## shape’s final geometry
            print("Line: ", x, "length: ", geom.length())
        else:
            x = geom.asMultiPolyline() ## shape’s final geometry
            print("MultiLine: ", x, "length: ", geom.length()) # print x 可能會資料量太大
    elif geom.type() == QgsWkbTypes.PolygonGeometry:
        if geomSingleType:
            x = geom.asPolygon()## shape’s final geometry
            print("Polygon: ", x, "Area: ", geom.area())
        else:
#            x = geom.asMultiPolygon()## shape’s final geometry
            print("Area: ", geom.area())
#            print("MultiPolygon: ", x, "Area: ", geom.area()) # 印了當機
    else:
        print("Unknown or invalid geometry")

    attrs = feature.attributes() # attributes
    print(attrs)
    break
除了使用feature['NAME']，也可以使用feature[index]，e.g. feature[1]。

Selection

uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

########################################################
import random as rd
hexcolor = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd', 'e', 'f']
def randomColor1():
    """
    * 使用#rrggbb形式
    """
    co = "#"
    for i in range(6):
        co += rd.choice(hexcolor)
    return co

def randomColor2():
    """
    * 使用(R, G, B, Alpha)形式
    """
    return QColor(rd.randint(0,255),rd.randint(0,255),rd.randint(0,255),rd.randint(0,255))
    

# 改變selection color
#iface.mapCanvas().setSelectionColor(QColor(250, 150, 100, 50)) # last number indicates alpha (transparency 0~255), optional
#iface.mapCanvas().setSelectionColor(QColor("#abcdef"))
#iface.mapCanvas().setSelectionColor(QColor("blue"))
#iface.mapCanvas().setSelectionColor(QColor(randomColor1()))
iface.mapCanvas().setSelectionColor(randomColor2())

ocean.selectAll() ## 選擇全部 see also: Select Features (or Edit > Select > Select Feature(s))

countries.selectByExpression('"fid"=\'186\' or "fid"=\'154\'', QgsVectorLayer.SetSelection) ## 根據條件選擇

addFeatures = [1, 2, 8, 9, 20, 25, 30, 100]
countries.select(addFeatures)

ocean.removeSelection() ## To remove the selection

iterating over selected features

歷遍selected feature。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

########################################################

countries.selectByExpression('"fid"=\'186\' or "fid"=\'154\'', QgsVectorLayer.SetSelection) ## 根據條件選擇

## 使用selectedFeatures()來取得某layer的selecction
selection = countries.selectedFeatures()
for feature in selection:
    print(feature['NAME'])

iterating over a subset of features

取得某區域內的元件(feature)。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

########################################################
areaOfInterest = QgsRectangle(-9.5,-32.2, 40, 36) # 選擇左下角及右上角座標
request = QgsFeatureRequest().setFilterRect(areaOfInterest) # 過濾該矩形內容
## 使用getFeatures(request)來取得request內的feature
for feature in countries.getFeatures(request):
    print(feature['NAME'])

Styling

利用PyQGIS來改變feature的style。首先選擇一個點的layer開啟: 
ulayer = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_50m_geography_regions_points", "twcity", "ogr")
ulayer.renderer().symbol().setSize(5)
ulayer.renderer().symbol().setColor(QColor("red"))
ulayer.renderer().symbol().symbolLayer(0).setShape(QgsSimpleMarkerSymbolLayerBase.Star)
ulayer.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(ulayer.id())
renderer()讓我們得到renderer物件(e.g. Single Symbol (QgsSingleSymbolRenderer), Categorized (QgsCategorizedSymbolRenderer), Graduated (QgsGraduatedSymbolRenderer)等)

Renderer

當我們提交(render)一個圖層，其外觀由圖層的renderer跟symbols給定。Symbols classes用來負責繪出元件，而renderers則決定那一個symbol用在哪一個feature。 
cities = QgsVectorLayer("C:/Maps/twcity.shp", "cities", "ogr") ## remember to use / insteas of \
if not cities.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(cities)
################################################################################
cities_renderer = cities.renderer() ## 取得layer cities的renderer
print(cities_renderer.type()) # > singleSymbol

## renderer type 有許多，可以使用以下方式列出:
print(QgsApplication.rendererRegistry().renderersList())

## 使用dump()方法來取得rederer內容
print(cities_renderer.dump()) # > SINGLE: MARKER SYMBOL (1 layers) color 114,155,111,255

## 先使用createSimple()函數來建立點(QgsMarkerSymbol)、線(QgsLineSymbol)或面(QgsFillSymbol)的simbol
## 然後使用setSymbol()函數來改變symbol

symbol = QgsMarkerSymbol.createSimple({'name': 'square', 'color': 'red'})
symbol.setSize(10) 
cities.renderer().setSymbol(symbol)
# show the change
cities.triggerRepaint() ##重畫

設定layer feature, selected feature, and stroke顏色

uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)

########################################################
import random as rd
def randomColor2():
    """
    * 使用(R, G, B, Alpha)形式
    """
    return QColor(rd.randint(0,255),rd.randint(0,255),rd.randint(0,255),rd.randint(0,255))
color = randomColor2()
print(f"{color.name()}")

countries.selectByExpression('"fid"=\'186\' or "fid"=\'154\'', QgsVectorLayer.SetSelection) ## 根據條件選擇

## 使用selectedFeatures()來取得某layer的selecction
selection = countries.selectedFeatures()
countries.renderer().symbol().setColor(QColor("blue")) ## 設定layer中feature的顏色
iface.mapCanvas().setSelectionColor( QColor("red") ) ## 設定layer中selected feature的顏色
for feature in selection:
    # 取得highlight物件
    highlight = QgsHighlight(iface.mapCanvas(), feature, countries)
    # 設定feature的邊框顏色
    highlight.setColor(color)

修改圖層(Modifying Vector Layers)

首先先檢查有哪些操作可以使用。 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_admin_0_countries"
countries = iface.addVectorLayer(uri, "countries", "ogr")
KH = QgsVectorLayer("C:\Maps\KHRoads.shp", "KH", "ogr")
ocean = iface.addVectorLayer("C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_10m_ocean", "Ocean", "ogr")
if not KH.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(KH)
########################################################

caps_string = KH.dataProvider().capabilitiesString()
print(caps_string) # 印出所有可進行操作

## 確認某一操作是否有支援
caps = KH.dataProvider().capabilities()
if caps & QgsVectorDataProvider.DeleteFeatures:
    print("可以刪除元件")

建立新的vector圖層

我們可以使用layer > Create Layer > New Shapefile Layer...來建立新的圖層，使用PyQGIS的方法: 
# 建立點的圖層，名稱為layername，使用memory data provider
layerl = QgsVectorLayer("Point", "layername", "memory")

# 增加資料欄位
from qgis.PyQt.QtCore import QVariant
pr = layerl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String),
                  QgsField("age",  QVariant.Int),
                  QgsField("size", QVariant.Double)])
layerl.updateFields() # 更新欄位

# 加入一個點
f = QgsFeature()
f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(120.3,22.6))) # 點位置
f.setAttributes(["Ada L.", 2, 0.3]) # 設定該點屬性
pr.addFeature(f) # 將該點加入到layer的feature
layerl.updateExtents() # 更新內容
QgsProject.instance().addMapLayer(layerl) # 將layer加入到canvas

# 統計屬性

print("欄位(field)數:", len(pr.fields()))
print("元件(feature)數:", pr.featureCount())
e = layerl.extent()
print("Extent:", e.xMinimum(), e.yMinimum(), e.xMaximum(), e.yMaximum()) # 座標的極值
 
for f in layerl.getFeatures():
    print("Feature:", f.id(), f.attributes(), f.geometry().asPoint()) # 所有feature的屬性

# 若想增加屬性欄位

layerl.startEditing() # 開始編輯
 
new_field_name = 'New field' # 新欄位名稱
layerl.addAttribute(QgsField(new_field_name, QVariant.String)) # 加入屬性，此欄位資料為string
layerl.updateFields() # 更新欄位
 
for f in layerl.getFeatures():
    print("Feature:", f.id(), f.attributes(), f.geometry().asPoint())

# 更新資料
new_field_value = 'New value'
for f in layerl.getFeatures():
    f[new_field_name] = new_field_value # 每一比資料的新欄位資料都是new_field_value
    layerl.updateFeature(f) # 增加後更新
 
layerl.commitChanges() # 確認改變
 
for f in layerl.getFeatures():
    print("Feature:", f.id(), f.attributes(), f.geometry().asPoint())

iface.vectorLayerTools().stopEditing(layerl) # 停止編輯

## 修該資料的替代方法:
my_field_name = 'New field 2'
my_field_value = 'New value 2'
 
with edit(layerl): # 使用with，無須手動呼叫startEditing(), commitChanges(), stopEditing()等
    layerl.addAttribute(QgsField(my_field_name, QVariant.String))
    layerl.updateFields()
    for f in layerl.getFeatures():
        f[my_field_name] = my_field_value
        layerl.updateFeature(f)

如果增加的是線的圖層: 
# 建立線的圖層，名稱為layername，使用memory data provider
layerl = QgsVectorLayer("LineString", "layername", "memory")

# 增加資料欄位
from qgis.PyQt.QtCore import QVariant
pr = layerl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String),
                  QgsField("age",  QVariant.Int),
                  QgsField("length", QVariant.Double)])
layerl.updateFields() # 更新欄位

lon1 = 120.3
lat1 = 22.6
lon2 = 120.5
lat2 = 22.8
linePoints = [QgsPoint(lon1,lat1), QgsPoint(lon2,lat2)]

line = QgsGeometry.fromPolyline(linePoints)
f = QgsFeature()
f.setGeometry(line) # 加入線
f.setAttributes(["River View", 2, 33.3]) # 設定該線屬性
pr.addFeature(f) # 將該線加入到layer的feature
layerl.updateExtents() # 更新內容
QgsProject.instance().addMapLayer(layerl) # 將layer加入到canvas

with edit(layerl): # 使用with，無須手動呼叫startEditing(), commitChanges(), stopEditing()等
    layerl.addAttribute(QgsField("from_lon", QVariant.Double))
    layerl.addAttribute(QgsField("from_lat", QVariant.Double))
    layerl.addAttribute(QgsField("to_lon", QVariant.Double))
    layerl.addAttribute(QgsField("to_lat", QVariant.Double))
    layerl.updateFields()
    for fea in layerl.getFeatures():
        fea["from_lon"] = lon1
        fea["from_lat"] = lat1
        fea["to_lon"] = lon2
        fea["to_lat"] = lat2
        layerl.updateFeature(fea)

練習:嘗試在線的圖層加入多條線。 
import random as rd
# 建立線的圖層，名稱為layername，使用memory data provider
layerl = QgsVectorLayer("LineString", "layername", "memory")

# 增加資料欄位
from qgis.PyQt.QtCore import QVariant
pr = layerl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String),
                  QgsField("age",  QVariant.Int),
                  QgsField("length", QVariant.Double)])
layerl.updateFields() # 更新欄位

# 加入多條線
lon1 = 120.3
lat1 = 22.6
lon2 = 120.5
lat2 = 22.8
lon = [lon1, lon2]
lat = [lat1, lat2]

with edit(layerl): # 使用with，無須手動呼叫startEditing(), commitChanges(), stopEditing()等
    layerl.addAttribute(QgsField("from_lon", QVariant.Double))
    layerl.addAttribute(QgsField("from_lat", QVariant.Double))
    layerl.addAttribute(QgsField("to_lon", QVariant.Double))
    layerl.addAttribute(QgsField("to_lat", QVariant.Double))
    
for i in range(3):
    linePoints = [QgsPoint(lon1,lat1), QgsPoint(lon2,lat2)]

    lon1, lon2 = lon2, lon2 + rd.random()*(1)
    lat1, lat2 = lat2, lat2 + rd.random()*(1)
    
    lon.append(lon2)
    lat.append(lat2)
    
    line = QgsGeometry.fromPolyline(linePoints)
    f = QgsFeature()
    f.setGeometry(line) # 加入線
    f.setAttributes(["River View", 2, 33.3, lon1, lat1, lon2, lat2]) # 設定該線屬性
    pr.addFeature(f) # 將該線加入到layer的feature
    layerl.updateExtents() # 更新內容
    QgsProject.instance().addMapLayer(layerl) # 將layer加入到canvas

練習:嘗試建立面的圖層。

Geometry

Geometry就是讓我們操作幾何圖形，原則上就是指點線面。在QGIS中，使用QgsGeometry class來儲存與表達。有時候一個geometry實際上是一堆單一的geometries的集合(Collection)，此類稱之為multi-part geometry，如果是包含單一類型，我們可以稱之為 multi-point, multi-linestring or multi-polygon。e.g. 一個區域由多個島嶼組成，則為multi-polygon。
在QGIS中，有多種建立geometry的方式，首先先回顧一下下例: 
# 建立點的圖層，名稱為layername，使用memory data provider
layer = QgsVectorLayer("Point", "layername", "memory")

# 增加資料欄位
from qgis.PyQt.QtCore import QVariant
pr = layer.dataProvider()
pr.addAttributes([QgsField("fid", QVariant.Int),
                  QgsField("longitude",  QVariant.Double),
                  QgsField("latitude", QVariant.Double)])
layer.updateFields() # 更新欄位

# 加入一個點
f = QgsFeature()
f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(120.3,22.6))) # 點位置
f.setAttributes([0, 120.3, 22.6]) # 設定該點屬性
pr.addFeature(f) # 將該點加入到layer的feature
layer.updateExtents() # 更新內容
QgsProject.instance().addMapLayer(layer) # 將layer加入到canvas

管理圖層

管理圖層的幾個方法:

Map Canvas

畫布(Canvas)顯示包含圖層(layers)的地圖(map)，我們可以使用map tools(panning, zooming, identifying layers, measuring, vector editing and others)來操控。畫布由QgsMapCanvas class in the qgis.gui module實現，每次我們操作地圖工具(例如移動或縮放)，地圖會重新提交(render)成為一個圖片(使用QgsMapRendererJob物件)，然後顯示在畫布上。此外canvas items是用來額外顯示或強調內容。例如rubber bands (used for measuring, vector editing etc.) or vertex markers。

Expression & Filtering

使用函數

我們可以開啟GeoPackage並使用OGR來取得所有的檔案名稱，並依此判斷某圖層是否存在於此GeoPackage，並依此作為開啟圖層或顯示錯誤訊息之依據。 
from osgeo import ogr
layerDir = 'C:/Maps/packages/natural_earth_vector.gpkg' # natural_earth_vector.gpkg的路徑
layers = [lay.GetName() for lay in ogr.Open(layerDir)] # 路境內所有的檔案名
print(layers)
# 判斷某圖層是否存在
print("Whatever" in layers) # False
print('ne_10m_parks_and_protected_lands_area' in layers) # True

#################################################
## 依以上方式來判斷某圖層存在後再進行開啟，否則顯示錯誤訊息

layerName1 = "Whatever"
layerName2 = 'ne_110m_ocean'
if layerName1 in layers:
    iface.addVectorLayer(layerDir + "|layername=" + layerName1, layerName1, 'ogr')
else:
    print(f"Error: There is no layer named {layerName1} in {layerDir}")

if layerName2 in layers:
    iface.addVectorLayer(layerDir + "|layername=" + layerName2, layerName2, 'ogr')
else:
    print(f"Error: There is no layer named {layerName2} in {layerDir}")
將上述流程寫成一個函數供使用: 
from osgeo import ogr
layerDir = 'C:/Maps/packages/natural_earth_vector.gpkg' # natural_earth_vector.gpkg的路徑

def toAddLayer(layerdir, layername):
    layers = [lay.GetName() for lay in ogr.Open(layerdir)] # 路境內所有的檔案名
    if layername in layers:
        iface.addVectorLayer(layerDir + "|layername=" + layername, layername, 'ogr')
    else:
        print(f"Error: There is no layer named {layername} in {layerdir}")

toAddLayer(layerDir, "Whatever") # False
toAddLayer(layerDir, "ne_110m_land") # True
若是想要可以一次開啟多個圖層: 
from osgeo import ogr
layerDir = 'C:/Maps/packages/natural_earth_vector.gpkg' # natural_earth_vector.gpkg的路徑

def toAddLayer(layerdir:str, layername:str):
    """
    * 開啟一個圖層
    """
    layers = [lay.GetName() for lay in ogr.Open(layerdir)] # 路境內所有的檔案名
    if layername in layers:
        iface.addVectorLayer(layerDir + "|layername=" + layername, layername, 'ogr')
    else:
        print(f"Error: There is no layer named {layername} in {layerdir}")

def toAddMultiLayers(layerdir:str, names:list):
    """
    * 開啟多個圖層
    """
    for name in names:
        toAddLayer(layerdir, name) # 呼叫 toAddLayer

toAddMultiLayers(layerDir, ["Whatever", "ne_110m_land", 'ne_110m_ocean']) # False, True, True

Map Rendering and Layouts

Communicating with the user

與使用者互動的方式。

Network analysis

Network analysis可以用來建立網路並求解網路問題(目前只有Dijkstra’s algorithm)。所以主要使用步驟可分為

Processing Tools

使用Processing Tools來協助。首先看以下的code: 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_110m_populated_places"
result = processing.run("native:buffer", 
    {'INPUT':uri,'DISTANCE':10,'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,
    'DISSOLVE':False,'OUTPUT':'memory:'})
QgsProject.instance().addMapLayer(result['OUTPUT'])

buffer方法定義: This algorithm computes a buffer area for all the features in an input layer, using a fixed or dynamic distance.。 建立對話視窗: 
uri = "C:/Maps/packages/natural_earth_vector.gpkg|layername=ne_110m_populated_places"
result = processing.createAlgorithmDialog("native:buffer", 
    {'INPUT':uri,'DISTANCE':10,'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,
    'DISSOLVE':False,'OUTPUT':'memory:'})
result.show()
processing種類繁多，可以參考Processing providers and algorithms。

練習看看最短路徑的做法: 
uri = "C:\Maps\KHRoads.shp"
result = processing.run("qgis:shortestpathpointtopoint", 
    {'INPUT':uri, 'START_POINT':"120.28999,22.82482", 'END_POINT':"120.34166, 22.82857", 'OUTPUT': "memory:"})
QgsProject.instance().addMapLayer(result['OUTPUT'])
也可以使用對話視窗: 
uri = "C:\Maps\KHRoads.shp"
iface.addVectorLayer(uri, "khroads", "ogr")
result = processing.createAlgorithmDialog("qgis:shortestpathpointtopoint", 
    {'INPUT':uri, 'START_POINT':"120.28999,22.82482", 'END_POINT':"120.34166, 22.82857", 'OUTPUT': "memory:"})
result.show()
加入開啟地圖，不然地圖無法顯示。

Open Street Map

是開放軟體，可以使用Leaflet來操控。我們可以下載Leaflet來使用，也可以直接使用在CDN的Hosted Version，在download頁面可以找到最新版。將其複製，然後貼到html檔案的head內。如下: 
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.8.0/dist/leaflet.css" integrity="sha512-hoalWLoI8r4UszCkZ5kL8vayOGVae1oxXe/2A4AO6J9+580uKHDO3JdHb7NzwwzK5xr/Fs0W40kiNHxM9vyTtQ==" crossorigin="" />
<script src="https://unpkg.com/leaflet@1.8.0/dist/leaflet.js" integrity="sha512-BB3hKbKWOc9Ez/TAwyWxNXeoV9c1v6FIeYiBieIWkpLjauysF18NzgR1MBNBXf8/KABdlkX68nAhlwcDFLGPCQ==" crossorigin=""></script>
接下來在body內建立一個放置地圖的地方，如下: 
<div id="mapid"></div>
加上相關的CSS讓其顯示。 
body {
    margin: 0;
    padding: 0;
}
#mapid {
    width: 100vw;
    height: 100vh;
}
接下來加上如下的javascripe code即可。 
let map = L.map('mapid', {
    center: [22.757296287028296, 120.3376949843434],
    zoom: 17
});

L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors'
}).addTo(map);
然後在瀏覽器應該就可以看到地圖。 以上是基本的起手式，接著我們就可以來操作此地圖了。

Marker

我們可以在地圖上加上標記，使用的語法如下:

Layer

之前介紹的Marker與Popup屬於UI Layers，而之前顯示的地圖屬於raster layer。若想顯示不同型態外觀的地圖，可至Leaflet providers網站選取想要的style。或參考官網之Basemap Providers。
Vector Layer則包含點線面等元件。

Event

也可以進行事件處理，例如: 
map.on('click', (e)=>{console.log(e.latlng)});
注意若是點擊在前述有popup的marker上，會出現popup而不會得到座標。不過倒是可以使用popup來顯示座標，如下: 
map.on('click', (e)=>{
    L.popup().setLatLng(e.latlng).setContent(e.latlng.toString()).openOn(map);
});

Routing

欲在地圖上繪製最短路徑路線並計算路線長度，可以使用套件(Plugin)，此處介紹Leaflet Routing Machine。做法如下，首先在head加上leaflet-routing-machine.css與leaflet-routing-machine.js: 
<link rel="stylesheet" href="https://unpkg.com/leaflet-routing-machine@latest/dist/leaflet-routing-machine.css" />
<script src="https://unpkg.com/leaflet-routing-machine@latest/dist/leaflet-routing-machine.js"></script>
若想要下載亦可，方式詳見官網。接下來在javascript內加上以下程式碼即可。 
L.Routing.control({
waypoints: [
    L.latLng(22.757296287028296, 120.3376949843434),
    L.latLng(22.795367, 120.337632)
]
}).addTo(map);
只要給兩地的經緯度即可計算並顯示。做個routing的例子: 
let map = L.map('mapid').setView([22.757296287028296, 120.3376949843434], 12);
L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors'
}).addTo(map);


L.Routing.control({
    waypoints: [
        L.latLng(22.757296287028296, 120.3376949843434),
        L.latLng(22.755367, 120.297632),
        L.latLng(22.727087, 120.310142),
        L.latLng(22.744107, 120.335462),
        L.latLng(22.757296287028296, 120.3376949843434)
    ]
}).addTo(map);
再試一個可以自訂marker與popup的例子: 
let map = L.map('mapid').setView([22.757296287028296, 120.3376949843434], 12);
L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors'
}).addTo(map);

L.Routing.control({
    waypoints: [
        L.latLng(22.66863, 120.30982),
        L.latLng(22.64067, 120.308404),
        L.latLng(22.757296287028296, 120.3376949843434)
    ],
    routeWhileDragging: true,
    lineOptions: { // 修改線型
        styles: [{ color: 'green', opacity: 1, weight: 5 }]
    },

    createMarker: function (i, waypoint, n) { // 自訂marker
        
        let marker = L.marker(waypoint.latLng, {
            draggable: true, // 可拖曳
            bounceOnAdd: false,
            bounceOnAddOptions: {
                duration: 1000,
                height: 800,
                function() {
                    (bindPopup(popup).openOn(map))
                }
            },
            icon: L.icon({
                iconUrl: 'airplane.png',
                iconSize: [30, 30],
                iconAnchor: [15, 15],
                popupAnchor: [0, 0],
                shadowUrl: '',
                shadowSize: [0, 0],
                shadowAnchor: [0, 0]
            })
            
        });
        //建立popup
        let popup = L.popup().setContent(`<span style="font-weight: bolder; color: blue;">${i}</span>:${waypoint.latLng}->${n}`);
        
        return marker.bindPopup(popup);
    }
}).addTo(map);