Add a field which is holding the area of each record (the workflow)

This commit is contained in:
Alireza Adli 2024-03-12 16:20:55 -04:00
parent 6c6d46255c
commit ce515b5427
2 changed files with 56 additions and 28 deletions

View File

@ -4,10 +4,11 @@ Buildings Geometries.
Developed by Alireza Adli
"""
from qgis.core import QgsApplication, QgsVectorLayer, QgsProject, \
QgsProcessingFeedback
from qgis.core import QgsApplication, QgsField, QgsProject, \
QgsProcessingFeedback, QgsVectorLayer, QgsVectorDataProvider, \
QgsExpressionContext, QgsExpressionContextUtils
from qgis.analysis import QgsNativeAlgorithms
from qgis.PyQt.QtCore import QVariant
import processing
import os
from services_scripts.basic_functions import create_folders, \
@ -17,7 +18,7 @@ from services_scripts.basic_functions import create_folders, \
# This function loads a layer.
# It's used to count the data records
def load_layer(path, layer_name):
the_layer = QgsVectorLayer(path, layer_name, "ogr")
the_layer = QgsVectorLayer(path, layer_name, 'ogr')
if not the_layer.isValid():
print(f'{layer_name} failed to load!')
else:
@ -74,7 +75,7 @@ params_clipping_nrcan = {'INPUT': fixed_nrcan_layer,
'FILTER_EXPRESSION': '',
'FILTER_EXTENT': None,
'OUTPUT': clipped_nrcan_layer}
processing.run("native:clip", params_clipping_nrcan)
processing.run('native:clip', params_clipping_nrcan)
print(f'Clipping of {fixed_nrcan_name} is completed.')
# After each step, data records are counted
@ -90,7 +91,7 @@ print(f'{clipped_nrcan_name} data count: {clipped_nrcan_length}')
# (smoothly and efficiently)
params_create_index_nrcan = {'INPUT': clipped_nrcan_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_nrcan)
processing.run('native:createspatialindex', params_create_index_nrcan)
print(f'Creating spatial index for {clipped_nrcan_name} is completed.')
# Reading Shared platform of geospatial data
@ -102,7 +103,7 @@ geoindex, geoindex_name = load_layer(geoindex_layer, 'GeoIndex')
print(f'{geoindex_name} data count: {geoindex.featureCount()}')
params_create_index_geo = {'INPUT': geoindex_layer, 'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_geo)
processing.run('native:createspatialindex', params_create_index_geo)
# Fixing the GeoIndex layer geometries
print(f'Fixing {geoindex_name} geometries')
@ -112,11 +113,11 @@ fixed_geoindex = 'C:/Users/a_adli/PycharmProjects/' \
params_fixing_geoindex = {'INPUT': geoindex,
'METHOD': 0,
'OUTPUT': fixed_geoindex}
processing.run("native:fixgeometries", params_fixing_geoindex)
processing.run('native:fixgeometries', params_fixing_geoindex)
params_create_index_fixed_geo = {'INPUT': fixed_geoindex,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_fixed_geo)
processing.run('native:createspatialindex', params_create_index_fixed_geo)
fixed_geoindex_read, fixed_geoindex_name = \
load_layer(fixed_geoindex, 'Fixed NRCan')
@ -134,7 +135,7 @@ params = {'INPUT': fixed_geoindex,
'FILTER_EXPRESSION': '',
'FILTER_EXTENT': None,
'OUTPUT': clipped_geoindex_layer}
processing.run("native:clip", params)
processing.run('native:clip', params)
print(f'Clipping {fixed_geoindex_name} is completed.')
clipped_geoindex, clipped_geoindex_name = \
@ -144,7 +145,7 @@ print(f'{clipped_geoindex_name} data count: {clipped_geoindex.featureCount()}')
# Creating spatial index fo the GeoIndex layer
params_create_index_clipped_geo = {'INPUT': clipped_geoindex_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_clipped_geo)
processing.run('native:createspatialindex', params_create_index_clipped_geo)
print(f'Creating spatial index for {clipped_geoindex_name} is completed.')
# Reading the Property Assessment (AKA uniteevaluationfonciere) dataset
@ -161,7 +162,7 @@ print(f'{property_assessment_layer_name} '
# Creating spatial index for the GeoIndex layer
params_create_index_property_assessment = {'INPUT': property_assessment_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex",
processing.run('native:createspatialindex',
params_create_index_property_assessment)
print(f'Creating spatial index for '
f'{property_assessment_layer_name} is completed.')
@ -194,12 +195,12 @@ for each in range(num_layers):
f'{(each + 1) * splitting_intervals}\r\n',
'OUTPUT': output_layer_path}
processing.run("native:extractbyexpression", splitting_params)
processing.run('native:extractbyexpression', splitting_params)
# Creating spatial index for each new layer
params_create_index_layers = {'INPUT': output_layer_path,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_layers)
processing.run('native:createspatialindex', params_create_index_layers)
# Putting the remaining features in one last layer
# (So, overall we're going to have [num_layers + 1] layers).
@ -217,12 +218,12 @@ last_splitting_params = {'INPUT': clipped_nrcan_layer,
f'$id >= {num_layers * splitting_intervals}\r\n',
'OUTPUT': output_layer_path}
processing.run("native:extractbyexpression", last_splitting_params)
processing.run('native:extractbyexpression', last_splitting_params)
# Create spatial index for the last partition (layer)
params_create_index_last = {'INPUT': output_layer_path,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_last)
processing.run('native:createspatialindex', params_create_index_last)
# Now we clip the property assessment with each of the partitions (layers)
# from the last process. The process outputs the pairwise clipped version of
@ -253,14 +254,14 @@ for each in range(num_clipped_layers):
'FILTER_EXPRESSION': '',
'FILTER_EXTENT': None,
'OUTPUT': clipped_layer}
processing.run("native:clip", pairwise_clipping_params)
processing.run('native:clip', pairwise_clipping_params)
# Creating spatial index for each new layer
params_create_index_clips = {'INPUT': clipped_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_clips)
processing.run('native:createspatialindex', params_create_index_clips)
print("Pairwise Clipping is completed.")
print('Pairwise Clipping is completed.')
# Finally, we merge all the clipped layers to make the
# pairwise clipped Property Assessment layer
@ -277,12 +278,12 @@ merging_params = {'LAYERS': pairwise_clipped_layers_path,
'CRS': None,
'OUTPUT': pairwise_clipped_property_assessment_layer}
processing.run("native:mergevectorlayers", merging_params)
processing.run('native:mergevectorlayers', merging_params)
# Create spatial index for the pairwise clipped Property Assessment
params_create_index_pairwised = \
{'INPUT': pairwise_clipped_property_assessment_layer, 'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_pairwised)
processing.run('native:createspatialindex', params_create_index_pairwised)
pairwise_clipped_property_assessment, \
pairwise_clipped_property_assessment_name = \
@ -319,7 +320,7 @@ processing.run('native:joinattributesbylocation',
params_create_index_joined_nrcan = \
{'INPUT': property_assessment_nrcan_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", params_create_index_joined_nrcan)
processing.run('native:createspatialindex', params_create_index_joined_nrcan)
property_assessment_nrcan, property_assessment_nrcan_name = \
load_layer(property_assessment_nrcan_layer,
@ -350,7 +351,7 @@ processing.run('native:joinattributesbylocation',
create_index_params = \
{'INPUT': property_assessment_nrcan_geo_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", create_index_params)
processing.run('native:createspatialindex', create_index_params)
property_assessment_nrcan_geo, property_assessment_nrcan_geo_name = \
load_layer(property_assessment_nrcan_geo_layer,
@ -374,11 +375,11 @@ property_nrcan_geo_deleted_duplicates_layer = \
'output_data/deleted_duplicates_property_and_all/delete_duplicates.shp'
delete_duplicates_params = {'INPUT': property_assessment_nrcan_geo_layer,
'OUTPUT': property_nrcan_geo_deleted_duplicates_layer}
processing.run("native:deleteduplicategeometries", delete_duplicates_params)
processing.run('native:deleteduplicategeometries', delete_duplicates_params)
create_index_params = {'INPUT': property_nrcan_geo_deleted_duplicates_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", create_index_params)
processing.run('native:createspatialindex', create_index_params)
property_nrcan_geo_deleted_duplicates, \
property_nrcan_geo_deleted_duplicates_name = \
@ -392,17 +393,44 @@ singled_parts_layer = \
'hydroquebec_archetype_gispy/data/' \
'output_data/singled_parts/singled_parts.shp'
# Multipart to single parts process: enables clicking for even small geometries
# that used to be, mistakenly, only clickable with a bigger geometry together
singled_parts_params = {'INPUT': property_nrcan_geo_deleted_duplicates_layer,
'OUTPUT': singled_parts_layer}
processing.run("native:multiparttosingleparts", params)
processing.run('native:multiparttosingleparts', params)
create_index_params = {'INPUT': singled_parts_layer,
'OUTPUT': 'Output'}
processing.run("native:createspatialindex", create_index_params)
processing.run('native:createspatialindex', create_index_params)
singled_parts, singled_parts_layer_name = \
load_layer(singled_parts_layer,
'Property Assessment Layer Joined with NRCan')
print(f'{singled_parts_layer_name} data count: {singled_parts.featureCount()}')
# Add an 'Area' field to the layer, to hold the area of each record
functionalities = singled_parts.dataProvider().capabilities()
if functionalities & QgsVectorDataProvider.AddAttributes:
new_field = QgsField('Area', QVariant.Double)
singled_parts.dataProvider().addAttributes([new_field])
singled_parts.updateFields()
# Retrieve the area of each record and put it in the new added field: 'Area'
singled_parts.startEditing()
idx = singled_parts.fields().indexFromName('Area')
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(singled_parts))
for feature in singled_parts.getFeatures():
area = feature.geometry().area()
feature[idx] = area
singled_parts.updateFeature(feature)
singled_parts.commitChanges()
qgs.exitQgis()

View File

@ -9,6 +9,6 @@ layer = QgsVectorLayer(
functionalities = layer.dataProvider().capabilities()
if functionalities & QgsVectorDataProvider.AddAttributes:
new_field = QgsField("Area", QVariant.Double)
new_field = QgsField('Area', QVariant.Double)
layer.dataProvider().addAttributes([new_field])
layer.updateFields()