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)

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.

Where do you generate the random distance and direction? Or are you expecting someone else to code this for you? The way things work here, you make a best effort and we point out the mistake, but if you don't make any effort, why should we? — Vince, Nov 02 '23 at 13:41

Taras · Accepted Answer · 2023-11-23T21:02:40.270

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.

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, 2pi)
    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 and get the output that will look like this:

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

References:

Thank you so much - I am really grateful. I run the script on my dataset (see dataset link in post) but I get the error "AttributeError: 'list' object has no attribute 'crs'". But the dataset has crs. See error output here https://mapsafe.xyz/masking_output.png — Pankaj Sharma, Nov 02 '23 at 23:23
you probably forgot [0] in input_layer = QgsProject.instance().mapLayersByName("all_clusters_kamloops")[0] — Taras, Nov 03 '23 at 08:38

Adding simple Geomasking (Donut masking) feature using PyQGIS

1 Answers1

Solution 1 : Projected Coordinate System

It uses translate() method of the QgsGeometry class

referring to the original point layer

creating a new layer for the output

accessing output layer provider

inheriting data from the original layer

processing new features

adding new layer to the map