GEOG 489
Advanced Python Programming for GIS

4.13.1 Lesson 4 Practice Exercise Solution and Explanation

PrintPrint

Solution

import qgis 
import qgis.core
import sys, os 
import pandas as pd 


# read data into pandas data frame 
data = pd.read_csv(r"C:\489\L4\exercise\L4exercise_data.csv") 

# class definition for PointObject 
class PointObject(): 

    def __init__(self, lat, lon):  
        self.lon = lon 
        self.lat = lat 

    def toQgsFeature(self):  
         feat = qgis.core.QgsFeature() 
         feat.setGeometry(qgis.core.QgsGeometry.fromPointXY(qgis.core.QgsPointXY(self.lon, self.lat))) 
         return feat 

firstObservations = []          # for storing objects of class pointWithID 
firstObservationsFeatures = []  # for storing objects of class QgsFeature 

# code for creating QgsApplication and initializing QGIS environment 
qgis_prefix = os.getenv("QGIS_PREFIX_PATH")                
qgis.core.QgsApplication.setPrefixPath(qgis_prefix, True) 
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis() 
 
# definition of class PointWithID derived from PointObject  
# to represent animal observation from the input data 
class PointWithID(PointObject): 
 
    def __init__(self, pID, lat, lon): 
        super(PointWithID, self).__init__(lat, lon) 
        self.pID = pID  # instance variable for storing animal ID 

    # overwriting the == operator to be based on the animal ID 
    def __eq__(self, other): 
        return self.pID == other.pID 

    # overwriting this method to include animal ID as attribute of QgsFeature created 
    def toQgsFeature(self): 
         feat = super(PointWithID, self).toQgsFeature() 
         feat.setAttributes([self.pID]) 
         return feat 

# create list of PointWithID object with first observations for each animal in the data frame 
for row in data.itertuples(index=False): 
    pointWithID = PointWithID(row[0], row[1], row[2])  
    if not pointWithID in firstObservations: # here __eq__() is used to do the comparison 
        firstObservations.append(pointWithID) 

# list comprehension for creating list of features from firstObservations list 
firstObservationsFeatures = [ o.toQgsFeature() for o in firstObservations ] 

# create new point layer with field for animal ID 
layer = qgis.core.QgsVectorLayer("Point?crs=EPSG:4326&field=AnimalID:string(255)", 'animal first observations' ,'memory') 

# add features to new layer 
prov = layer.dataProvider() 
prov.addFeatures(firstObservationsFeatures) 

# save layer as GeoPackage file 
qgis.core.QgsVectorFileWriter.writeAsVectorFormat( layer, r"C:\489\L4\exercise\firstobservations.gpkg", "utf-8",layer.crs(), "GPKG") 

# clean up 
qgs.exitQgis() 

Explanation

  • In line 34 the new class PointWithID is declared to be derived from class PointObject.
  • The constructor in line 36 calls the constructor of the base class and then adds the additional instance variable for the animal ID (called pID).
  • __eq__() is defined in lines 43 and 44 to compare the pID instance variables of both involved objects and return True if they are equal.
  • In lines 48 and 49, we call toQgsFeature() of the base class PointObject and then take the returned feature and add the animal ID as the only attribute.
  • In lines 53 to 56, we use itertuples() of the data frame to loop through its rows. We create a new PointWithID object in variable pointWithID from the cell values of the row and then test whether an object with that ID is already contained in list firstObservations. If not, we add this object to the list.
  • In line 59, we use toQgsFeature() to create a QgsFeature object from each PointWithID object in the list. These are then added to the new layer in line 66.
  • When creating the new layer in line 62, we have to make sure to include the field for storing the animal ID in the string for the first parameter.