Skip to content

Commit

Permalink
Merge pull request #287 from osundwajeff/points_per_grid_cell
Browse files Browse the repository at this point in the history 
Points per grid cell
  • Loading branch information
@osundwajeff
osundwajeff authored Sep 20, 2024
2 parents 4a03d11 + 8e5b1b7 commit e9ca92b
Show file tree
Hide file tree
Showing 4 changed files with 325 additions and 87 deletions.
184 changes: 98 additions & 86 deletions src/qgis_gender_indicator_tool/jobs/create_grids.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,91 +31,103 @@ def create_grids(self, layer, output_dir, crs, merged_output_path):
# Check if the merged grid already exists
if os.path.exists(merged_output_path):
print(f"Merged grid already exists: {merged_output_path}")
return QgsVectorLayer(
merged_output_path, "merged_grid", "ogr"
) # Load the existing merged grid layer

layer = QgsVectorLayer(layer, "country_layer", "ogr")
if not layer.isValid():
raise ValueError("Invalid country layer")

# Reproject the country layer if necessary
if layer.crs() != crs:
layer = processing.run(
"native:reprojectlayer",
{
"INPUT": layer,
"TARGET_CRS": crs,
"OUTPUT": "memory:",
},
feedback=QgsProcessingFeedback(),
)["OUTPUT"]

all_grids = []

# Loop through each feature in the polygon layer
for feature in layer.getFeatures():
geom = feature.geometry()

# Check if the geometry is multipart
if geom.isMultipart():
parts = (
geom.asGeometryCollection()
) # Separate multipart geometry into parts
else:
parts = [geom] # Single part geometry

# Loop through each part of the geometry
for part_id, part in enumerate(parts):
part_area = part.area()

# Get the extent of each part
part_extent = part.boundingBox()

# Define the output grid path for each part
grid_output_path = (
f"{output_dir}/grid_{feature.id()}_part_{part_id}.gpkg"
)

# Check if the grid already exists
if os.path.exists(grid_output_path):
print(f"Grid file already exists: {grid_output_path}")
grid_layer = QgsVectorLayer(
grid_output_path, "grid_layer", "ogr"
) # Load the existing grid layer
return merged_output_path
else:
layer = QgsVectorLayer(layer, "country_layer", "ogr")
if not layer.isValid():
raise ValueError("Invalid country layer")

# Reproject the country layer if necessary
if layer.crs() != crs:
layer = processing.run(
"native:reprojectlayer",
{
"INPUT": layer,
"TARGET_CRS": crs,
"OUTPUT": "memory:",
},
feedback=QgsProcessingFeedback(),
)["OUTPUT"]

all_grids = []

# Loop through each feature in the polygon layer
for feature in layer.getFeatures():
geom = feature.geometry()

# Check if the geometry is multipart
if geom.isMultipart():
parts = (
geom.asGeometryCollection()
) # Separate multipart geometry into parts
else:
print(f"Creating grid: {grid_output_path}")
# Define grid creation parameters
grid_params = {
"TYPE": 2, # Rectangle (polygon)
"EXTENT": part_extent, # Use the extent of the current part
"HSPACING": self.h_spacing, # Horizontal spacing
"VSPACING": self.v_spacing, # Vertical spacing
"CRS": crs, # Coordinate reference system (CRS)
"OUTPUT": grid_output_path, # Output path for the grid file
}

# Create the grid using QGIS processing
grid_result = processing.run("native:creategrid", grid_params)
grid_layer = grid_result["OUTPUT"] # Get the grid layer

# Clip the grid to the polygon feature (to restrict it to the boundaries)
clipped_grid_output_path = (
f"{output_dir}/clipped_grid_{feature.id()}_part_{part_id}.gpkg"
parts = [geom] # Single part geometry

# Loop through each part of the geometry
for part_id, part in enumerate(parts):
part_area = part.area()

# Get the extent of each part
part_extent = part.boundingBox()

# Define the output grid path for each part
grid_output_path = (
f"{output_dir}/grid_{feature.id()}_part_{part_id}.shp"
)
clip_params = {
"INPUT": grid_layer, # The grid we just created
"OVERLAY": layer, # The layer we're clipping to
"OUTPUT": clipped_grid_output_path,
}
clip_result = processing.run("native:clip", clip_params)
grid_layer = clip_result["OUTPUT"] # The clipped grid

# Add the generated or loaded grid to the list
all_grids.append(grid_layer)

# Merge all grids into a single layer
print(f"Merging grids into: {merged_output_path}")
merge_params = {"LAYERS": all_grids, "CRS": crs, "OUTPUT": merged_output_path}
merged_grid = processing.run("native:mergevectorlayers", merge_params)["OUTPUT"]
return merged_grid

# Check if the grid already exists
if os.path.exists(grid_output_path):
print(f"Grid file already exists: {grid_output_path}")
grid_layer = QgsVectorLayer(
grid_output_path, "grid_layer", "ogr"
) # Load the existing grid layer
# Clip the grid to the polygon feature (to restrict it to the boundaries)
clipped_grid_output_path = f"{output_dir}/clipped_grid_{feature.id()}_part_{part_id}.shp"
clip_params = {
"INPUT": grid_layer, # The grid we just created
"OVERLAY": layer, # The layer we're clipping to
"OUTPUT": clipped_grid_output_path,
}
clip_result = processing.run("native:clip", clip_params)
grid_layer = clip_result["OUTPUT"] # The clipped grid
else:
print(f"Creating grid: {grid_output_path}")
# Define grid creation parameters
grid_params = {
"TYPE": 2, # Rectangle (polygon)
"EXTENT": part_extent, # Use the extent of the current part
"HSPACING": self.h_spacing, # Horizontal spacing
"VSPACING": self.v_spacing, # Vertical spacing
"CRS": crs, # Coordinate reference system (CRS)
"OUTPUT": grid_output_path, # Output path for the grid file
}

# Create the grid using QGIS processing
grid_result = processing.run("native:creategrid", grid_params)
grid_layer = grid_result["OUTPUT"] # Get the grid layer

# Clip the grid to the polygon feature (to restrict it to the boundaries)
clipped_grid_output_path = f"{output_dir}/clipped_grid_{feature.id()}_part_{part_id}.shp"
clip_params = {
"INPUT": grid_layer, # The grid we just created
"OVERLAY": layer, # The layer we're clipping to
"OUTPUT": clipped_grid_output_path,
}
clip_result = processing.run("native:clip", clip_params)
grid_layer = clip_result["OUTPUT"] # The clipped grid

# Add the generated or loaded grid to the list
all_grids.append(grid_layer)

# Merge all grids into a single layer
print(f"Merging grids into: {merged_output_path}")
merge_params = {
"LAYERS": all_grids,
"CRS": crs,
"OUTPUT": merged_output_path,
}
merged_grid = processing.run("native:mergevectorlayers", merge_params)[
"OUTPUT"
]

return merged_grid
139 changes: 139 additions & 0 deletions src/qgis_gender_indicator_tool/jobs/points_per_grid_cell.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
import os
from qgis.PyQt.QtCore import QVariant
from qgis.core import (
QgsVectorLayer,
QgsField,
QgsSpatialIndex,
QgsProcessingFeedback,
)
import processing
from .create_grids import GridCreator
from .extents import Extents


class RasterPointGridScore:
def __init__(self, country_boundary, pixel_size, output_path, crs, input_points):
self.country_boundary = country_boundary
self.pixel_size = pixel_size
self.output_path = output_path
self.crs = crs
self.input_points = input_points

def raster_point_grid_score(self):
"""
Generates a raster based on the number of input points within each grid cell.
:param country_boundary: Layer defining the country boundary to clip the grid.
:param cellsize: The size of each grid cell.
:param output_path: Path to save the output raster.
:param crs: The CRS in which the grid and raster will be projected.
:param input_points: Layer of point features to count within each grid cell.
"""

# Create grid
self.h_spacing = 100
self.v_spacing = 100
create_grid = GridCreator(h_spacing=self.h_spacing, v_spacing=self.v_spacing)
output_dir = os.path.join("output")
merged_output_path = os.path.join(
output_dir, "merged_grid.shp"
) # Use Shapefile

# Create grid layer using Shapefile
grid_layer = create_grid.create_grids(
self.country_boundary, output_dir, self.crs, merged_output_path
)
grid_layer = QgsVectorLayer(merged_output_path, "merged_grid", "ogr")

# Add score field
provider = grid_layer.dataProvider()
field_name = "score"
if not grid_layer.fields().indexFromName(field_name) >= 0:
provider.addAttributes([QgsField(field_name, QVariant.Int)])
grid_layer.updateFields()

# Create spatial index for the input points
# Reproject the country layer if necessary
if self.input_points.crs() != self.crs:
self.input_points = processing.run(
"native:reprojectlayer",
{
"INPUT": self.input_points,
"TARGET_CRS": self.crs,
"OUTPUT": "memory:",
},
feedback=QgsProcessingFeedback(),
)["OUTPUT"]
point_index = QgsSpatialIndex(self.input_points.getFeatures())

# Count points within each grid cell and assign a score
reclass_vals = {}
for grid_feat in grid_layer.getFeatures():
grid_geom = grid_feat.geometry()
# Get intersecting points
intersecting_points = point_index.intersects(grid_geom.boundingBox())
num_points = len(intersecting_points)

# Reclassification logic: assign score based on the number of points
if num_points >= 2:
reclass_val = 5
elif num_points == 1:
reclass_val = 3
else:
reclass_val = 0

reclass_vals[grid_feat.id()] = reclass_val

# Apply the score values to the grid
grid_layer.startEditing()
for grid_feat in grid_layer.getFeatures():
grid_layer.changeAttributeValue(
grid_feat.id(),
provider.fieldNameIndex(field_name),
reclass_vals[grid_feat.id()],
)
grid_layer.commitChanges()

merged_output_vector = os.path.join(
output_dir, "merged_grid_vector.shp"
) # Use Shapefile for merged output

# Merge grids into a single Shapefile layer
Merge = processing.run(
"native:mergevectorlayers",
{"LAYERS": [grid_layer], "CRS": None, "OUTPUT": "memory:"},
)

merge = Merge["OUTPUT"]

extents_processor = Extents(
output_dir, self.country_boundary, self.pixel_size, self.crs
)

# Get the extent of the vector layer
country_extent = extents_processor.get_country_extent()
xmin, ymin, xmax, ymax = (
country_extent.xMinimum(),
country_extent.yMinimum(),
country_extent.xMaximum(),
country_extent.yMaximum(),
)

# Rasterize the clipped grid layer to generate the raster
rasterize_params = {
"INPUT": merge,
"FIELD": field_name,
"BURN": 0,
"USE_Z": False,
"UNITS": 1,
"WIDTH": self.pixel_size,
"HEIGHT": self.pixel_size,
"EXTENT": f"{xmin},{xmax},{ymin},{ymax}",
"NODATA": -9999,
"OPTIONS": "",
"DATA_TYPE": 5, # Use Int32 for scores
"OUTPUT": self.output_path,
}

processing.run(
"gdal:rasterize", rasterize_params, feedback=QgsProcessingFeedback()
)
2 changes: 1 addition & 1 deletion test/test_create_grids.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ def test_create_grids(self):
os.path.dirname(__file__), "test_data/admin/Admin0.shp"
)
self.output_dir = os.path.join(os.path.dirname(__file__), "output")
self.merged_output_path = os.path.join(self.output_dir, "merged_grid.gpkg")
self.merged_output_path = os.path.join(self.output_dir, "merged_grid.shp")
self.utm_crs = QgsCoordinateReferenceSystem("EPSG:32620") # UTM Zone 20N
self.h_spacing = 100
self.v_spacing = 100
Expand Down
Loading

0 comments on commit e9ca92b

Please sign in to comment.