Rivers as polygones and treated as barriers do not end up with continuous stretches of NoData

[nicolasray]

Some river data sets are polygons rather than lines, and usually the ones in polygones are of much higher precision and are used more frequently. However, when rasterized in AccessMod to be treated as barriers, they can translate in a non-continuous stretch of NoData (which is not the case when lines are rasterized), which permits passage when it should not! See below an example of it, showing the merged landcover where the polygone rivers (in blue) were treated as barriers:

The solution could be in the settings of the rasterized function used in AccessMod: if all the cells touching the polygones were considered as the river this could solve the issue.

A current manual fix could be to get the centerline of the polygon rivers as an additional line shapefile to add as barrier during the landcover merged. This way we are ensured there is no passage possible. A function to do that in ArcGIS Pro is here: https://pro.arcgis.com/en/pro-app/tool-reference/topographic-production/polygon-to-centerline.htm

I have not fund a working solution in QGIS (the QGIS plugin HCMGIS that is supposed to do it is extremely slow and I couldn't get a result from it).

[fxi]

First thought : why not use a r.mask on polygon barrier, fill inner cell as 0 or 1 and continue with the script ? Not sur if r.mask behave like v.to.rast. The flag -d works for lines and boundary, so we should rasterize twice, once as polygon, once with densified line on boundary version of the polygon.

Refs : r.mask v.to.rast

[SteeveEbener]

I think it would be important to first identify the cause before trying to apply a remedy: why is this happening with the polygon and not with the center line?

There is maybe an issue with the rasterisation of polygons in AccessMod.

Applying the method proposed by Nicolas (converting all the cells touching the polygon as NoData) could significantly expand the water bodies area depending on the resolution of the input data => this is more a path than a real solution.

[fxi]

So, there is no generic solution to this issue.

As long as we use raster based barriers, there will be issues with thin/densified lines. There is no right answer. Rasterizing at low resolution produces artefacts.

An unsupervised polygon centerline algorithm could cover some cases, but I'm not sure that this approach could fit all cases. This should probably be done by the user outside AccessMod.

Additionally – for water bodies at least – I'm pretty sure that there is already a lot of data available with the centerline already computed. It's not uncommon data. Then, the user can always add as many polygons and lines she wants to build the merged land cover layer : there is no limit on the number of stack items used.

However. As we have already agreed to use densified lines for barriers of type 'lines', I've added an option to produce the same effect for polygons. Any cell touching polygons will be used as a barrier. There will be islands produced from polygon with this method. If no facilities are mistakenly placed in the middle of those areas, the travel time algorithm will never visit those locations. For the cases where islands are not wanted, a second stack item could simply be generated to fill those islands.

Here is an example of barrier using the densified polygon edges approach :

This change is available now in testing (92b01c8afb)

[SteeveEbener]

Is it me or the figure included by Fred in his last message shows that only the border of the polygons (lines) have been rasterized and not its whole surface?

Is it possible to get an answer to the question I raised in my previous message and to have a figure showing the result without the densified polygon edges approach before this get implemented in AccessMod?

Thanks

[fxi]

@SteeveEbener : I wrote a paragraph that describes the strategy along the figure you mentioned, I invite you to read it. In few words, the new option produces impassable barriers for the travel time algorithm but also islands (see below, fig. 4). Remaining islands can be removed easily with an additional step, aka "two passes rasterisation" (see below, fig. 5).

Then, the answer to your question is also written in my comment, but maybe you are not familiar with basic rasterisation process.

If you want to learn more about the tool I use in the background, please read the link I've provided v.to.rast

Here are all possible cases I've mentioned :

Line rasterisation using thin line strategy fig.1 : One adjacent cell selected along all segments. Produce bridges.

Line rasterisation using densified line: fig.2 : All touching cells selected along all segments [ updated see comment by @nicolasray ]

Polygon rasterisation (current behaviour) fig.3 : All cells where the center point is under a polygon are selected. Produce bridges. see @nicolasray comment.

Polygon edges rasterisation using densified line (new option) fig.4 : All touching cells selected along all segments. [ updated see comment by @nicolasray ]

Polygon rasterisation + polygon edges rasterisation using densified line ( need 2 passes ) fig.5 : All touching cells selected along all segments + All cells where the center point is under a polygon are selected [ updated see comment by @nicolasray ]

Is that more clear to you ?

[nicolasray]

Something is confusing in the previous post: Figure 2 (Line rasterization using densified line) is NOT the current behaviour in AccessMod when one rasterizes lines. Looking at several "landcover merged" example I have, this is rather what happen with the current code:

Fred, can you confirm? And this is the correct way (and the default way of rasterizing lines in many GIS package). It should stay like that and we should not use the densified method for lines (unless the densified method is indeed what is currently implemented, which I think it is..).

Now for polygones, it is quite clear for me that the correct method should be the one in Figure 5 (if the rasterization of the line segments along the polygone is the same as the one now implemented in AccessMod). Deciding that all cells in which a river line (or a river/lake bank in case of polygone) passes through should be treated as barrier makes sense. We could decide otherwise (e.g. some land in a cell is more important than some water), but how to justify it? I'm really in favor of adopting the behaviour of Figure 5, as it will avoid many potential errors. But I'm also in favor of having this behaviour used as an option (unchecked by default), with some good explanation and example in the user guide as to when to use it.

[fxi]

@nicolasray Oups, you are right, sorry. I've updated the figure.

For polygons, I've implemented fig.4, but I can patch resulting grid of fig.3 to fig.4, that would gives fig.5. It's easy. Or I can implement both.

[fxi]

@nicolasray and yes, I confirm, densified line is the only way to avoid artificial bridges, and this is the approach used in AccessMod 5. As already said, this is a workaround to the fact that we don't use a system that would process vector barriers in r.walk.accessmod – the module that produces travel time analysis.

[nicolasray]

OK, great. Can you also update Figure 4 and 5: there is one pixel on the second row from top that is not touching the line. I will use your Figures in the manual later when we explain this situation. But let's see what Steeve has to say first.

[fxi]

@nicolasray I've moved the line and forgot to update the cells. Thanks for noticing this.

If you want the source OpenOffice document (spreadsheet), here it is : issue_259_rasterisation_drawing.ods.zip You can change the colors: it's simply conditional styling.

[SteeveEbener]

Hi Fred,

I have been doing rasterisation for a long time but thanks for the refresher. It would indeed be good to have these figures in the user manual.

I read your paragraph. The figure you used was just not illustrating what you would expect as the final result for the rasterisation of a polygon and the reason for my additional questions.

The solution is therefore indeed to implement the approach resulting in figure 5 here above, method which is similar to the one included after figure 5 (coincidence...) in section 3.3.1.1 of the user manual.

=> it is great to have this directly implemented in AccessMod. Thank you very much for that.

[fxi]

@SteeveEbener : Yes, fig.5 is indeed the method described in the section 3.3.1.1 of the manual.

@nicolasray I've found a way to extract polygons skeleton :

fig. 6 : Extracting skeleton of polygons using voronoi diagram approach

I tried to use the method described in the manual of v.voronoi. The polygons have been properly dissolved first.

Should I implement this instead ?

[nicolasray]

Yes, we should implement this centerline method instead, as it ensures that we have a continuous stretch of data, and ensure we don't have the extra NoData cells that were appearing in the method "polygon edges rasterisation using densified line".

[fxi]

So, I've implemented this.

I've also added a step to create a buffer of 0 meter to dissolve quickly all polygon into one. So, the user can input any polygon layer without worrying about bad/unexpected center lines in nested polygons.

That's great, but there is edges cases with the 'center line' approach.

Using the default parameters (thin=-1), a single centerline is produced (fig. 6, left) but it would be more appropriate to use a skeleton (fig.6, right).

fig. 6 centerline vs skeleton

With real data, I've got this:

 _fig. 7 `thin = -1` -> pure center line_

 _fig. 8 `thin = map resolution` -> skeleton_

 _fig. 9 `thin = 2 * map resolution ` -> skeleton_

 _fig. 10 `thin = 10 * map resolution ` -> skeleton_

Should be let the user set its own value of thin ?

I will push 'centerline' aka 'thin=-1' in branch testing now and you can play with it. Then next week we can proceed to choose the best approach.

[nicolasray]

Good to see these alternatives. Of course, it's only relevant to use the rasterized centerline (likely the one Figure 9, but we need to discuss to understand what these alternatives mean in different cases. ) if it is a second pass on top of a first pass that needs to use the v.to.rast (so that we don't miss the many barrier pixels in the core areas of polygones). I'm assuming you have implemented these two passes (with the second one currently "thin = 1") in the testing branch.

[fxi]

thin is an option of v.voronoi. It's done before v.to.rast.

You can just look at the code here: amServer.R#L1430

The fig.9 is using a maximum dangle line length of 10 * project resolution. In that case – 10km. Then, it's rasterized using densified lines in v.to.rast as any other line. It's not a centerline.

The fig.7 is a centerline.

The different alternatives are giving more or less importance towards minor features' shape.

[SteeveEbener]

I think it would be important to see the result of the implementation of both "passes" on different size of polygons for a given resolution to see if this combination is indeed covering the whole area of the original polygon or if the polygon border as a line does also need to be included in the set before the rasterization.

[nicolasray]

After some discussions, it appears there would be no fit-all solution, so we decided to implement the following on any polygon shapefile used in the landcover merged:

• the current behaviour (using v.to.rast) is always applied
• a second pass using the "v.voronoi" is offered as an option to the user to apply. If chosen, the user must also specify the value of the parameter "thin" (see FIg. 7-10 in above comment), which will give either a pure centerline added, or some other "skeleton" that will give more accurate results in some cases. Clear explanations (with figures) will be added in the user guide to explain this.

[SteeveEbener]

Thanks for looking further into this issue.

If we go this way, which I agree is a good one, we should then also give the option for the boundaries of the polygons to be converted into lines (polylines) and then rasterized.

This is the option that we have been documenting, and therefore promoting, in the user manual since version 4.

This could be a separated option of this second pass and could be either applied on its own or in combination with the "v.voronoi".

[nicolasray]

We can indeed add a second option to additionally add the polygon boundaries converted as barrier pixels. However, note that this is not what has been promoted in the manual (it seems we never went to this details of explanation).

[SteeveEbener]

Section 3.3.1.1 of version 5.0 user manual (Section 3.2.1.1 for version 4.0):

"In addition to the above, the rasterization of elongated water body polygons (and other barriers to movement stored as polygons) does not always generate continuous surfaces, unlike when rasterizing vector lines. Figure 5 below provides an illustration of such problem. As we can see in Figure 5c, working with lower resolution (larger grid sizes) will create discontinuities, therefore generating artificial bridges in these areas.

The solution to address this problem is to:

1. Convert your polygon water body layer into polyline.
2. Merge the polyline resulting from step 1 with your line format barrier layer (see section 3.3.1.8 for details).
3. Use both the original polygon layer and the merged line layer as barrier when creating the combined land cover layer in AccessMod."

[fxi]

I wrote a script that seems to be quite robust to solve a lot of cases while remaining quite fast : 4.5 minutes for 77 million cells at 100 meter resolution. It's not perfect, but seems more reliable and less impacting than other solutions. It's available as an option in the testing branch, version 5.6.26

@YG1975 : You can do some testing with your project, if you want. It's not ready for production, though.

Description of the process

Here are the steps used. The sample is a thin area of a river in BFA.

1. Select a base resolution, lower than the one from the project, and set the computational grid based on that value. Example, 20 meters instead of 100. It could be 10 meters. The number of cells in the computational grid will be much larger.
2. Convert the polygon barrier to raster.
3. Add a buffer to the vectorized barrier large enough to cover neighborhood cells in the resolution of the project. If 100 m, a buffer of 100/2 m should do that.
4. Reset the grid resolution to the original resolution.
5. Extract linear components using a tool like r.thin. This will compute the skeleton of the areas.
6. Vectorize the output to obtain lines from the skeleton.
7. Rasterize the lines using the densification option.

Various cases

• Complex lake
• Close water bodies
• Very thin features
• Intricate river

Comparison with the 'polygon to lines' version

• Various comparison with the 'polygon to lines' version – described in the manual (yellow = rasterized polylines, red = r.thin result using the new script.)

[SteeveEbener]

Hi Fred,

Thanks for these different tests. This is very useful.

Just one point of clarification: we are anyway talking about combining the result of the new method (whatever it is) with the result of the rasterization of the polygon itself, right? It should be the case according to me.

If you are in agreement, would it be possible to include the cells resulting from the rasterization of the polygons in the screenshot comparing the two approaches?

Many thanks in advance for that and take care

[fxi]

For now, the new option create 2 layers : one barrier as lines (skeleton) and one barrier as areas (polygon rasterized). The user can inspect them in the raster preview or in a GIS independently, and, in the stack, merge them.

[fxi]

I missed the second part. I can include some cases with areas and skeletons combined. I can do it on Friday or the 10th of January.

[SteeveEbener]

Please also include the same cases with areas and "Polygon to polyline" combined to be able to make the comparison. Thanks

We should discuss if not combining the two layers by default is a good option because this could result in population located on real barriers and/or fake bridges if the user decide not to use both of them.

[nicolasray]

We have now made additional testing on the latest implementation of the new option. The results are shown below on a variety of polygon cases: the white areas represent the rasterisation of the barriers as it appears after the landcover merged tool, taking into account (in the merging process) one barrier layer as lines (skeleton) and one barrier layer as areas (polygon rasterized). Theses areas represent well the vector polygons (in red in figures below), without any fake bridges of cells. There are two parameters for this skeleton function: 1. Polygons rasterisation resolution, and 2. Polygon skeleton buffer size. How we set these two parameters could change a bit how the rasterisation is done. To allow a maximum of flexibility to the users, we give default values for these two parameters to give good results in most cases, ut the user can change these values if necessary. In the user interface, they will appear as follow:

This is now implemented in the "testing" version and will be pushed in the main branch very soon.

Test images are below:

[SteeveEbener]

Thanks for the update.

Can we see the screenshot presenting the results over the same areas without checking the "Add polygon skeleton as line barriers(s) for comparison? Such comparison will also be useful for the user manual.

I also think that the text below the two options should be adjusted for the user to better understand what this means. It is for example not clear to me that the polygon rasterization resolution parameter is part of the skeleton process

[nicolasray]

We have adjusted the text in the skeleton module. For the requested screenshots, we'll do it as soon as possible.

[SteeveEbener]

I just looked at the skeleton component of the barrier menu and while the text for the "Polygons rasterisation resolution [m]" field has been modified, the one for the "Polygon skeleton buffer size [m]" one remains the same => could we please also include a short intro sentence that describes what this field is about? Thanks