{"id":13036,"date":"2022-08-30T07:23:00","date_gmt":"2022-08-30T05:23:00","guid":{"rendered":"https:\/\/www.opengis.ch\/?p=13036"},"modified":"2022-10-13T09:07:04","modified_gmt":"2022-10-13T07:07:04","slug":"writing-a-feature-based-processing-algorithm-at-the-example-of-m-value-interpolation","status":"publish","type":"post","link":"https:\/\/www.opengis.ch\/de\/2022\/08\/30\/writing-a-feature-based-processing-algorithm-at-the-example-of-m-value-interpolation\/","title":{"rendered":"Writing a feature-based processing algorithm at the example of M-value interpolation"},"content":{"rendered":"\n

<\/p>\n\n\n\n

Amongst all the processing algorithms already available in QGIS, sometimes the one thing you need is missing. <\/p>\n\n\n\n

This happened not a long time ago, when we were asked to find a way to continuously visualise traffic on the Swiss motorway network<\/strong> (polylines) using frequently measured traffic volumes from discrete measurement stations (points) alongside the motorways. In order to keep working with the existing polylines, and be able to attribute more than one value of traffic to each feature, we chose to work with the M-values. M-values are a per-vertex attribute like X, Y or Z coordinates. They contain a measure value, which typically represents time or distance. But they can hold any numeric value.<\/p>\n\n\n\n

In our example, traffic measurement values are provided on a separate point layer and should be attributed to the M-value of the nearest vertex of the motorway polylines. Of course, the motorway features should be of type LineStringM<\/em> in order to hold an M-value. We then should interpolate the M-values for each feature over all vertices<\/strong> in order to get continuous values along the line (i.e. a value on every vertex). This last part is not yet existing as a processing algorithm in QGIS.<\/p>\n\n\n\n

This article describes how to write a feature-based processing algorithm<\/strong> based on the example of M-value interpolation along LineStrings. <\/p>\n\n\n\n

Feature-based processing algorithm<\/h2>\n\n\n\n

The pyqgis class QgsProcessingFeatureBasedAlgorithm<\/code><\/a> <\/strong>is described as follows: \u201cAn abstract QgsProcessingAlgorithm base class for processing algorithms which operates \u201cfeature-by-feature\u201d. \u00a0<\/em><\/p>\n\n\n\n

Feature based algorithms are algorithms which operate on individual features in isolation. These are algorithms where one feature is output for each input feature, and the output feature result for each input feature is not dependent on any other features present in the source. [\u2026]<\/em><\/p>\n\n\n\n

Using QgsProcessingFeatureBasedAlgorithm as the base class for feature based algorithms allows shortcutting much of the common algorithm code for handling iterating over sources and pushing features to output sinks. It also allows the algorithm execution to be optimised in future (for instance allowing automatic multi-thread processing of the algorithm, or use of the algorithm in \u201cchains\u201d, avoiding the need for temporary outputs in multi-step models).<\/em>\u201d<\/p>\n\n\n\n

In other words, when connecting several processing algorithms one after the other – e.g. with the graphical modeller – these feature-based processing algorithms can easily be used to fill in the missing bits. <\/p>\n\n\n\n

Compared to the standard QgsProcessingAlgorithm<\/code><\/a><\/strong> the feature-based class implicitly iterates over each feature when executing and avoids writing wordy loops explicitly fetching and applying the algorithm to each feature.\u00a0<\/p>\n\n\n\n

Just like for the QgsProcessingAlgorithm<\/code> (a template can be found in the Processing Toolbar > Scripts > Create New Script from Template<\/em>), there is quite some boilerplate code in the QgsProcessingFeatureBasedAlgorithm<\/code>. The first part is identical to any QgsProcessingAlgorithm<\/code>.<\/p>\n\n\n\n

After the description of the algorithm (name, group, short help, etc.), the algorithm is initialised with def initAlgorithm<\/code><\/em><\/strong>, defining input and output.\u00a0<\/p>\n\n\n\n

<\/p>\n\n\n\n

In our M-value example:<\/em><\/p>\n\n\n\n

    def initAlgorithm(self, config=None):\n        self.addParameter(\n            QgsProcessingParameterFeatureSource(\n                self.INPUT,\n                self.tr('Input layer'),\n                [QgsProcessing.TypeVectorAnyGeometry]\n            )\n        )\n        self.addParameter(\n            QgsProcessingParameterFeatureSink(\n                self.OUTPUT,\n                self.tr('Output layer')\n            )\n        )<\/code><\/pre>\n\n\n\n
While in a regular processing algorithm<\/strong> now follows def processAlgorithm(self, parameters, context, feedback)<\/code><\/strong>, in a feature-based algorithm<\/strong> we use def processFeature(self, feature, context, feedback)<\/code><\/em><\/strong>. This implies applying the code in this block to each feature of the input layer.\u00a0<\/p>\n\n\n\n
\n
! Do not use def processAlgorithm<\/code><\/strong><\/em> in the same script, otherwise your feature-based processing algorithm will not work !<\/td><\/tr><\/tbody><\/table><\/figure>\n<\/div><\/div>\n\n\n\n
Interpolating M-values<\/h2>\n\n\n\n
This actual processing part can be copied and added almost 1:1 from any other independent python script, there is little specific syntax to make it a processing algorithm. Only the first line below really. <\/p>\n\n\n\n
In our M-value example:<\/em><\/p>\n\n\n\n
    def processFeature(self, feature, context, feedback):\n        \n        try:\n            geom = feature.geometry()\n            line = geom.constGet()\n            vertex_iterator = QgsVertexIterator(line)\n            vertex_m = []\n\n            # Iterate over all vertices of the feature and extract M-value\n\n            while vertex_iterator.hasNext():\n                vertex = vertex_iterator.next()\n                vertex_m.append(vertex.m())\n\n            # Extract length of segments between vertices\n\n            vertices_indices = range(len(vertex_m))\n            length_segments = [sqrt(QgsPointXY(line[i]).sqrDist(QgsPointXY(line[j]))) \n                for i,j in itertools.combinations(vertices_indices, 2) \n                if (j - i) == 1]\n\n            # Get all non-zero M-value indices as an array, where interpolations \n              have to start\n\n            vertex_si = np.nonzero(vertex_m)[0]\n            \n            m_interpolated = np.copy(vertex_m)\n\n            # Interpolate between all non-zero M-values - take segment lengths between \n              vertices into account\n\n            for i in range(len(vertex_si)-1):\n                first_nonzero = vertex_m[vertex_si[i]]\n                next_nonzero = vertex_m[vertex_si[i+1]]\n                accum_dist = itertools.accumulate(length_segments[vertex_si[i]\n                                                                  :vertex_si[i+1]])\n                sum_seg = sum(length_segments[vertex_si[i]:vertex_si[i+1]])\n                interp_m = [round(((dist\/sum_seg)*(next_nonzero-first_nonzero)) + \n                            first_nonzero,0) for dist in accum_dist]\n                m_interpolated[vertex_si[i]:vertex_si[i+1]] = interp_m\n\n            # Copy feature geometry and set interpolated M-values, \n              attribute new geometry to feature\n\n            geom_new = QgsLineString(geom.constGet())\n            \n            for j in range(len(m_interpolated)):\n                geom_new.setMAt(j,m_interpolated[j])\n                \n            attrs = feature.attributes()\n            \n            feat_new = QgsFeature()\n            feat_new.setAttributes(attrs)\n            feat_new.setGeometry(geom_new)\n\n        except Exception:\n            s = traceback.format_exc()\n            feedback.pushInfo(s)\n            self.num_bad += 1\n            return []\n        \n        return [feat_new]\n<\/code><\/pre>\n\n\n\n
In our example, we get the feature\u2019s geometry, iterate over all its vertices (using the QgsVertexIterator<\/code><\/a>) and extract the M-values as an array. This allows us to assign interpolated values where we don’t have M-values available. Such missing values are initially set to a value of 0 (zero).<\/p>\n\n\n\n
We also extract the length of the segments between the vertices. By gathering the indices of the non-zero M-values of the array, we can then interpolate between all non-zero M-values<\/strong>, considering the length that separates the zero-value vertex from the first and the next non-zero vertex.<\/p>\n\n\n\n
For the iterations over the vertices to extract the length of the segments between them as well as for the actual interpolation between all non-zero M-value vertices we use the library itertools<\/a><\/em>. This library provides different iterator building blocks that come in quite handy for our use case. <\/p>\n\n\n\n
Finally, we create a new geometry by copying the one which is being processed and setting the M-values to the newly interpolated ones.<\/p>\n\n\n\n
And that’s all there is really!<\/p>\n\n\n\n
\n
Alternatively, the interpolation can be made using the interp<\/code> function of the numpy<\/code> library. Some parts where our manual<\/em> method gave no values, interp.numpy<\/code><\/a> seemed more capable of interpolating. It remains to be judged which version has the more realistic<\/em> results.<\/td><\/tr><\/tbody><\/table><\/figure>\n<\/div><\/div>\n\n\n\n
<\/p>\n\n\n\n
Styling the result via M-values<\/h2>\n\n\n\n
The last step is styling our output layer in QGIS, based on the M-values (our traffic M-values are categorised from 1 [a lot of traffic -> dark red] to 6 [no traffic -> light green]). This can be achieved by using a Single Symbo<\/em>l symbology with a Marker Line<\/em> type „on every vertex“<\/strong>. As a marker type, we use a simple round point. Stroke style<\/em> is „no pen“ and Stroke fill<\/em> is based on an expression:<\/p>\n\n\n\n
with_variable(\n\n'm_value', m(point_n($geometry, @geometry_point_num)),\n\n\tCASE WHEN @m_value = 6\n\t\tTHEN color_rgb(140, 255, 159)\n\n\t\tWHEN @m_value = 5\n\t\t\tTHEN color_rgb(244, 252, 0)\n\n\t\tWHEN @m_value = 4\n\t\t\tTHEN color_rgb(252, 176, 0)\n\n\t\tWHEN @m_value = 3\n\t\t\tTHEN color_rgb(252, 134, 0)\n\n\t\tWHEN @m_value = 2\n\t\t\tTHEN color_rgb(252, 29, 0)\n\n\t\tWHEN @m_value = 1\n\t\t\tTHEN color_rgb(140, 255, 159)\n\n\t\tELSE\n\t\t\tcolor_hsla(0,100,100,0)\n\n\tEND\n)<\/code><\/pre>\n\n\n\n
And voil\u00e0! Wherever we have enough measurements on one line feature, we get our motorway network continuously coloured according to the measured traffic volume.<\/p>\n\n\n\n
\n
\"\"
Motorway network – the different lanes are regrouped for each direction. M-values of the vertices closest to measurement points are attributed the measured traffic volume. The vertices are coloured accordingly.<\/figcaption><\/figure>\n\n\n\n
\"\"
Trafic on motorway network after „manual“ M-value interpolation. <\/figcaption><\/figure>\n\n\n\n
\"\"
Trafic on motorway network after M-value interpolation using numpy.<\/figcaption><\/figure>\n<\/figure>\n\n\n\n
<\/p>\n\n\n\n
One disclaimer at the end<\/em>: We get this seemingly continuous styling only because of the combination of our „complex“ polylines (containing many vertices) and the zoomed-out view of the motorway network. Because really, we’re styling many points and not directly the line itself. But in our case, this is working very well.<\/p>\n\n\n\n
If you’d like to make your custom processing algorithm available through the processing toolbox in your QGIS, just put your script in the folder containing the files related to your user profile:<\/p>\n\n\n\n
profiles > default > processing > scripts <\/code><\/pre>\n\n\n\n
You can directly access this folder by clicking on Settings<\/em> > User Profiles<\/em> > Open Active Profile Folder<\/em> in the QGIS menu. <\/p>\n\n\n\n
That way, it’s also available for integration in the graphical modeller.<\/p>\n\n\n\n
\"\"
Extract of the Graphical<\/em> Modeler<\/em> sequence. „Interpolate M-values neg“ refers to the custom feature-based processing algorithm described above.<\/p><\/figcaption><\/figure>\n\n\n\n

You can download the above-mentioned processing scripts (with numpy and without numpy) here<\/a><\/strong>.<\/p>\n\n\n\n
Happy processing!<\/p>\n","protected":false},"excerpt":{"rendered":"
Amongst all the processing algorithms already available in QGIS, sometimes the one thing you need is missing.  This happened not a long time ago, when we were asked to find a way to continuously visualise traffic on the Swiss motorway network (polylines) using frequently measured traffic volumes from discrete measurement […]<\/p>\n","protected":false},"author":38,"featured_media":13060,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_themeisle_gutenberg_block_has_review":false,"jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[36,14,15,161],"tags":[125],"class_list":["post-13036","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-processing","category-python","category-qgis","category-uncategorised","tag-qgis-org"],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2022\/08\/Screenshot-2022-07-26-at-10.25.35.png?fit=1592%2C1186&ssl=1","jetpack-related-posts":[{"id":13880,"url":"https:\/\/www.opengis.ch\/de\/2023\/09\/05\/analyzing-and-visualizing-large-scale-fire-events-using-qgis-processing-with-st-dbscan\/","url_meta":{"origin":13036,"position":0},"title":"Analyzing and visualizing large-scale fire events using QGIS processing with ST-DBSCAN","author":"Mathieu","date":"5. September 2023","format":false,"excerpt":"A while back, one of our ninjas added a new algorithm in QGIS\u2019 processing toolbox named ST-DBSCAN Clustering, short for spatio temporal density-based spatial clustering of applications with noise. The algorithm regroups features falling within a user-defined maximum distance and time duration values. This post will walk you through one\u2026","rel":"","context":"In "GIS"","block_context":{"text":"GIS","link":"https:\/\/www.opengis.ch\/de\/category\/gis\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2023\/09\/96.png?fit=755%2C566&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2023\/09\/96.png?fit=755%2C566&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2023\/09\/96.png?fit=755%2C566&ssl=1&resize=525%2C300 1.5x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2023\/09\/96.png?fit=755%2C566&ssl=1&resize=700%2C400 2x"},"classes":[]},{"id":3911,"url":"https:\/\/www.opengis.ch\/de\/2018\/05\/28\/how-to-filter-features-in-qgis-using-the-graphical-processing-modeler\/","url_meta":{"origin":13036,"position":1},"title":"How to filter features in QGIS using the graphical processing modeler","author":"Matthias Kuhn","date":"28. Mai 2018","format":false,"excerpt":"This article describes a new algorithm for the processing modeler called feature filter algorithm. If you are already familiar with ETL concepts and the graphical modeler, you can directly jump to the section the feature filter algorithm. Building workflows for repetitive tasks When building workflows for simple or complex geodata\u2026","rel":"","context":"In "Expressions"","block_context":{"text":"Expressions","link":"https:\/\/www.opengis.ch\/de\/category\/programming\/expressions\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/05\/modeller.png?fit=1016%2C497&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/05\/modeller.png?fit=1016%2C497&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/05\/modeller.png?fit=1016%2C497&ssl=1&resize=525%2C300 1.5x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/05\/modeller.png?fit=1016%2C497&ssl=1&resize=700%2C400 2x"},"classes":[]},{"id":3624,"url":"https:\/\/www.opengis.ch\/de\/2018\/06\/22\/threads-in-pyqgis3\/","url_meta":{"origin":13036,"position":2},"title":"Using Threads in PyQGIS3","author":"Marco Bernasocchi","date":"22. Juni 2018","format":false,"excerpt":"While porting a plugin to QGIS3 I decided to also move all it's threading infrastructure to QgsTasks. Here three possible variants to implement this.the first uses the static method QgsTask.fromFunction and is simpler to use. A great quick solution. If you want need control you can look at the second\u2026","rel":"","context":"In "PyQt"","block_context":{"text":"PyQt","link":"https:\/\/www.opengis.ch\/de\/category\/programming\/python\/pyqt\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/06\/pexels-wendy-van-zyl-1212179-scaled.jpg?fit=1200%2C800&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/06\/pexels-wendy-van-zyl-1212179-scaled.jpg?fit=1200%2C800&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/06\/pexels-wendy-van-zyl-1212179-scaled.jpg?fit=1200%2C800&ssl=1&resize=525%2C300 1.5x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/06\/pexels-wendy-van-zyl-1212179-scaled.jpg?fit=1200%2C800&ssl=1&resize=700%2C400 2x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/06\/pexels-wendy-van-zyl-1212179-scaled.jpg?fit=1200%2C800&ssl=1&resize=1050%2C600 3x"},"classes":[]},{"id":2235,"url":"https:\/\/www.opengis.ch\/de\/2017\/05\/02\/qgis-expressions-engine-performance-boost\/","url_meta":{"origin":13036,"position":3},"title":"QGIS Expressions Engine: Performance boost","author":"Matthias Kuhn","date":"2. Mai 2017","format":false,"excerpt":"Expressions in QGIS are more and more widely used for all kinds of purposes. For example the recently introduced geometry generators allow drawing awesome effects with modified feature geometries on the fly. The last days at the QGIS developer meeting 2017, I spent some time looking into and improving the\u2026","rel":"","context":"In "Expressions"","block_context":{"text":"Expressions","link":"https:\/\/www.opengis.ch\/de\/category\/programming\/expressions\/"},"img":{"alt_text":"","src":"","width":0,"height":0},"classes":[]},{"id":4330,"url":"https:\/\/www.opengis.ch\/de\/2018\/11\/06\/qml-widgets-qgis\/","url_meta":{"origin":13036,"position":4},"title":"The new QML widgets in QGIS – When widgets get unbridled","author":"Dave Signer","date":"6. November 2018","format":false,"excerpt":"Individuality is the definition of freedom. And freedom is the fundamental requirement of man\u2019s mind. QGIS possibly cannot give you all the freedom you require in life. But at least a lot of freedom in how you manage your work. QGIS 3.4.0 LTR was released last week and it comes\u2026","rel":"","context":"In "GIS"","block_context":{"text":"GIS","link":"https:\/\/www.opengis.ch\/de\/category\/gis\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/11\/pie_chart.png?fit=1067%2C521&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/11\/pie_chart.png?fit=1067%2C521&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/11\/pie_chart.png?fit=1067%2C521&ssl=1&resize=525%2C300 1.5x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/11\/pie_chart.png?fit=1067%2C521&ssl=1&resize=700%2C400 2x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2018\/11\/pie_chart.png?fit=1067%2C521&ssl=1&resize=1050%2C600 3x"},"classes":[]},{"id":13455,"url":"https:\/\/www.opengis.ch\/de\/2022\/11\/29\/qgis-relations-their-widgets-and-the-plugins-of-them\/","url_meta":{"origin":13036,"position":5},"title":"QGIS Relations, their Widgets and the Plugins of them","author":"Damiano","date":"29. November 2022","format":false,"excerpt":"This blog post is about QGIS relations and how they are edited in the attribute form with widgets in general, as well as some plugins that override the relations editor widget to improve usability and solve specific use cases. The start is quite basic. If you are already a relation\u2026","rel":"","context":"In "QGIS"","block_context":{"text":"QGIS","link":"https:\/\/www.opengis.ch\/de\/category\/gis\/qgis\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2022\/11\/feature.png?fit=641%2C463&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2022\/11\/feature.png?fit=641%2C463&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/www.opengis.ch\/wp-content\/uploads\/2022\/11\/feature.png?fit=641%2C463&ssl=1&resize=525%2C300 1.5x"},"classes":[]}],"jetpack_shortlink":"https:\/\/wp.me\/pbdBtI-3og","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/posts\/13036","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/users\/38"}],"replies":[{"embeddable":true,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/comments?post=13036"}],"version-history":[{"count":58,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/posts\/13036\/revisions"}],"predecessor-version":[{"id":13120,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/posts\/13036\/revisions\/13120"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/media\/13060"}],"wp:attachment":[{"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/media?parent=13036"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/categories?post=13036"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.opengis.ch\/de\/wp-json\/wp\/v2\/tags?post=13036"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}