PyQGIS – Adding Simple Geomasking (Donut Masking) Feature Using PyQGIS

affine transformationmaskingpointpyqgis

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.

Best Answer

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

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

References:

Related Solutions

Google Earth Engine – How to Mask Pixels Inside Polygons

There may be many ways to do that, but when I had to do it I did:

var maskInside = function(image, geometry) {
  var mask = ee.Image.constant(1).clip(geometry).mask().not()
  return image.updateMask(mask)
}

And put that function in a package (geetools for the code editor) that you can use in this way:

var tools = require('users/fitoprincipe/geetools:tools')  
var i = ee.Image.random(1)
var masked = tools.geometry.maskInside(i, geometry)
Map.addLayer(masked)

PyQGIS – Add LineString from Points

The fundamental reason why your code is not working is that you are trying to add a feature to your layer without making it editable. If you want to immediately commit the changes to the layer's underlying data store, you should work with layer.dataProvider() instead. Also, you should remove the last line from your code snippet, since the active layer is already added to the project, and replace it with a few lines to reflect the changes to the layer.

points = []
x = QgsPointXY(783648,1397243)
points.append(x)
y = QgsPointXY(788883,1396853)
points.append(y)

layer = iface.activeLayer()

if layer.dataProvider().capabilities() & QgsVectorDataProvider.AddFeatures:
    feature = QgsFeature(layer.fields())
    feature.setGeometry(QgsGeometry.fromPolylineXY(points))
    feature.setAttributes([0, 12])
    layer.dataProvider().addFeatures([feature])
    layer.updateExtents()
    layer.triggerRepaint()
    iface.mapCanvas().refresh()
    # Optional...
    iface.zoomToActiveLayer()

On the other hand, if you want the opportunity to undo/ redo then save the edits via the GUI, you can just call layer.startEditing() before adding the feature to the layer:

points = []
x = QgsPointXY(783648,1397243)
points.append(x)
y = QgsPointXY(788883,1396853)
points.append(y)

layer = iface.activeLayer()

if layer.type() == QgsMapLayerType.VectorLayer:
    feature = QgsFeature(layer.fields())
    feature.setGeometry(QgsGeometry.fromPolylineXY(points))
    feature.setAttributes([0, 12])

    if not layer.isEditable():
        layer.startEditing()
        
    layer.addFeature(feature)
    layer.updateExtents()
    layer.triggerRepaint()
    iface.mapCanvas().refresh()
    # Optional...
    iface.zoomToActiveLayer()

I strongly recommend reading the PyQGIS Cookbook section on Modifying Vector Layers.