Extract Z-values from multipoint feature in Pyqgis

Question

I have a program that stores elevation data in a multipoint vector using QgsPoint.

 caps3 = mp_layer3.dataProvider().capabilities() #check layer capabilities
 if caps3 & QgsVectorDataProvider.AddFeatures:
 for loop 1: (for loop over several lines)
 while loop 1: (while loop to save 3-D points on lines) 
 ... PROCESSING STUFF EXCLUDED FOR BREVITY ...
 mp.addGeometry(QgsPoint(pntxy.x(), pntxy.y(), zval[1], 0))
 
 xselevs.setGeometry(mp)
 xselevs.setAttributes([indx3])
 provider3.addFeature(xselevs)
 mp_layer3.updateExtents() 
 QgsProject.instance().addMapLayer(mp_layer3)

This all works great and I end up with a layer that has several multipoint features in it. I can access x, y, z, and m in QGIS without any problem. However, I want to be able to read back the layer in pyqgis, iterating over the features and reading the x, y, and z values. I came up with the following:

 XSPntFeats = mp_layer3.getFeatures() #get all cross sections in new XS point file
 for XSPnts in XSPntFeats: #loops through each multipoint contained in the shapefile
 for MPpart in XSPnts.geometry().asMultiPoint(): #loops through each part
 print("x-y-z: ", MPpart.x(), MPpart.y(), MPpart.z())

Works great for reading x and y, but alas, no z. It throws an error since z is not a part of the asMultiPoint() function. I can't find an iterator function that allows me to loop through each feature and read the points (including z) in them. I have done a lot of searching and fiddling and I cannot find a way to iterate over the features to get z back (or m for that matter).

Any ideas?

Answer 1

There are a few options for iterating through geometries, all fairly similar,a nd all of which will give you the same results - for MultiPoints. If you want to iterate through other geometry types, you will want to choose the right one for the job. For example, vertices() can be useful if you want to get the vertices of a MultiLineString or MultiPolygon without iterating through the parts.

1. Use QgsGeometry.vertices()

for f in vl.getFeatures():
 for v in f.geometry().vertices():
 print(v.x(), v.y(), v.z())

2. Use QgsGeometry.constGet()

for f in vl.getFeatures():
 for v in f.geometry().constGet():
 print(v.x(), v.y(), v.z())

3. Use QgsGeometry.constParts()

for f in vl.getFeatures():
 for v in f.geometry().constParts():
 print(v.x(), v.y(), v.z())

For the following raw data (ie a layer with two MultiPoints):

<QgsMultiPoint: MultiPointZ ((151.176021 -33.838767 100),(151.170861 -33.840853 101),(151.167229 -33.8419 102))>
<QgsMultiPoint: MultiPointZ ((151.140611 -33.807581 103),(151.156736 -33.792638 104))>

the output of each of the above is:

151.176021 -33.838767 100.0
151.170861 -33.840853 101.0
151.167229 -33.8419 102.0
151.140611 -33.807581 103.0
151.156736 -33.792638 104.0

To duplicate the data set up in QGIS, you can use the following PyQGIS code:

vl = QgsVectorLayer("MultiPointZ", "TempPoints", "memory")
mp1 = QgsMultiPoint()
mp2 = QgsMultiPoint()
mp1.addGeometry(QgsPoint(151.176021, -33.838767, 100))
mp1.addGeometry(QgsPoint(151.170861, -33.840853, 101))
mp1.addGeometry(QgsPoint(151.167229, -33.841900, 102))
mp2.addGeometry(QgsPoint(151.140611, -33.807581, 103))
mp2.addGeometry(QgsPoint(151.156736, -33.792638, 104))
vl.startEditing()
feat = QgsFeature()
feat.setGeometry(mp1)
vl.addFeature(feat)
feat.setGeometry(mp2)
vl.addFeature(feat)
vl.commitChanges()

Answer 2

I have three multipoints, first and second have three points each and third have only one.

How about:

for e, f in enumerate(QgsProject.instance().mapLayersByName("Draped")[0].getFeatures(),1): #For each feature
 print(f"Point: {e}")
 g = f.geometry() #Fetch its geometry
 for part in g.parts(): #For each part of the geometry
 print(round(part.z(),2)) #Each part is a QgsPoint Z, print the z value of it
#Point: 1
#684.04
#353.73
#470.05
#Point: 2
#603.41
#602.38
#552.46
#Point: 3
#311.85