Adding simple Geomasking (Donut masking) feature using PyQGIS

Question

I would like to add a simple halo masking feature/plugin in QGIS. Halo/donut masking basically transforms (shifts) all points in a layer to a random distance between two user-provided distances. During this process, these points are also transformed at random angles (auto-generated).

Once users have loaded the shapefile as a layer, they should be able to enter a minimum and maximum distance. From within this buffer, random distances should be generated for each point to be transformed. Also, a random angle should be generated for the transformation of each point. A new layer should be created containing the masked layer.

I have written the following code but need to develop it further. How can I transform the points into a new layer?

from PyQt4.QtCore import QVariant
from qgis.core import QGis, QgsFeature, QgsField, QgsFields
from processing.tools.vector import VectorWriter
import random
user_defined_minimum = 100
user_defined_maximum = 200
original_layer = iface.activeLayer()
for feature in original_layer.getFeatures():
 point = feature.geometry().asPoint()
 x, y = point.x(), point.y()
 print(str(x) + " " + str(y))
New_layer= processing.getObject(original_layer)
fields = QgsFields()
writer = processing.VectorWriter(output_file, None, fields.toList(), QGis.WKBMultiPolygon, New_layer.crs())
features = processing.features(New_layer)
for feat in features:
 x, y = random.uniform(1.5, 1.9) 
 writer.addFeature(feat)

enter image description here

This QGIS plugin has been requested since our MapSafe Data Sovereignty tool was released last month. MapSafe offers a complete approach for sovereign data owners to safeguard sensitive geospatial datasets by geomasking, encrypting, and notarising them from within the browser. The tool: https://www.mapsafe.xyz

More details are here: https://onlinelibrary.wiley.com/doi/epdf/10.1111/tgis.13094 (pdf)

The masking aspect is shown in a video here: https://www.mapsafe.xyz/safeguarding-guide.html Dataset link: https://mapsafe.xyz/all_clusters_kamloops.zip

In comparison to plugins on affine transformations, this plugin should be very simple.

Answer 1

Preamble: Unfortunatelly I did not read your articles, but I hope I understood your issue correctly.

Let's assume there is an input point layer called 'points2', see the image below.

input

The task is the same as yours to transform (shift) all points in a layer to a random distance between two user-provided distances and transform at a random angle (auto-generated).

Solution 1 : Projected Coordinate System

It uses translate() method of the QgsGeometry class

The code below also utilizes two Python packages: math and random.

In this example, I am working with EPSG:5348 which is a Projected Coordinate System.

Proceed with Plugins > Python Console > Show Editor and paste the script below:

# imports
from random import uniform
from math import pi, sin, cos
from qgis.core import QgsVectorLayer, QgsGeometry, QgsFeature, QgsProject
def geomasking(geom: QgsGeometry, min_distance: float = 0.0, max_distance: float = 10.0) -> QgsGeometry:
 """
 Only for Projected Coordinate System !
 Translates point geometry at a random angle and random distance within a range [min_distance : max_distance]
 Parameters:
 ==========
 :param geom: the feature's geometry
 :param min_distance: the minimum distance
 :param max_distance: the maximum distance
 """
 angle = uniform(0, 2*pi)
 distance = uniform(min_distance, max_distance)
 dx = distance * cos(angle)
 dy = distance * sin(angle)
 geom.translate(dx, dy)
 return geom
# referring to the original point layer
input_layer = QgsProject.instance().mapLayersByName("points2")[0]
# creating a new layer for the output
output_layer = QgsVectorLayer(f"Point?crs={input_layer.crs().authid()}&index=yes", "geomask", "memory")
# accessing output layer provider
provider = output_layer.dataProvider()
# inheriting data from the original layer
provider.addAttributes(input_layer.fields())
output_layer.updateFields()
# processing new features
for feat in input_layer.getFeatures():
 # providing geometry
 new_feat = QgsFeature()
 new_geom = geomasking(feat.geometry())
 new_feat.setGeometry(new_geom)
 # providing attributes
 new_feat.setAttributes(feat.attributes())
 # adding new feature to the output layer
 provider.addFeature(new_feat)
# adding new layer to the map
QgsProject.instance().addMapLayer(output_layer)

Solution 2 : Spherical Coordinate System

It uses computeSpheroidProject() method of the QgsDistanceArea class

The code below also utilizes two Python packages: math and random.

In this example, I am working with EPSG:4326 which is a Geographic/Spherical Coordinate System.

Proceed with Plugins > Python Console > Show Editor and paste the script below:

# imports
from math import pi
from random import uniform
from qgis.core import QgsVectorLayer, QgsGeometry, QgsFeature, QgsProject, QgsDistanceArea
def geomasking(geom: QgsGeometry, min_distance: float = 0.0, max_distance: float = 10.0) -> QgsGeometry:
 """
 Only for Spherical Coordinate System !
 Translates point geometry at a random angle and random distance within a range [min_distance : max_distance]
 Parameters:
 ==========
 :param geom: the feature's geometry
 :param min_distance: the minimum distance
 :param max_distance: the maximum distance
 """
 angle = uniform(0, 2 * pi)
 distance = uniform(min_distance, max_distance)
 da = QgsDistanceArea()
 da.setEllipsoid('WGS84')
 new_point = da.computeSpheroidProject(geom.asPoint(), distance, angle)
 return QgsGeometry.fromPointXY(new_point)
# referring to the original point layer
input_layer = QgsProject.instance().mapLayersByName("points2")[0]
# creating a new layer for the output
output_layer = QgsVectorLayer(f"Point?crs={input_layer.crs().authid()}&index=yes", "geomask", "memory")
# accessing output layer provider
provider = output_layer.dataProvider()
# inheriting data from the original layer
provider.addAttributes(input_layer.fields())
output_layer.updateFields()
# processing new features
for feat in input_layer.getFeatures():
 # providing geometry
 new_feat = QgsFeature()
 new_geom = geomasking(feat.geometry())
 new_feat.setGeometry(new_geom)
 # providing attributes
 new_feat.setAttributes(feat.attributes())
 # adding new feature to the output layer
 provider.addFeature(new_feat)
# adding new layer to the map
QgsProject.instance().addMapLayer(output_layer)

Change the inputs in the function if necessary (the default distance varies between 0 and 10). Press Run script run script and get the output that will look like this:

result

Addendum: Make yourself also familiar with Vincenty's formulae.

---

References:

• Rotating features with multiple geometries in one layer using PyQGIS
• Creating new layer with all attributes inherited from original layer using PyQGIS
• Azimuth calculations using PyQGIS
• Extracting x,y coordinates from QgsGeometry: Point(x, y) using PyQGIS
• Difference between QgsPoint, QgsPointXY and QgsGeometry.fromPointXY() in PyQGIS
• Vincenty or (Project + Euclidean) for accurate distances