Archives de catégorie : Billets

The power of array in qgis : match two sets of lines together

On Georezo, a French forum, someone asked how to retrieve the streets closest to pipelines.

The QGIS expression

So, let’s suppose we have one data for pipelines (source) and one for streets (target), here is a QGIS expression that retrieves the street data object index closest to each pipeline.

with_variable(
  'mean_distances',
  aggregate(
    'streets', 
    'array_agg', 
    array_avg(
      with_variable(
        'nodes',
        nodes_to_points(geometry(@parent)),
        array_foreach(
          generate_series(1, num_points(@nodes)),
          distance(point_n(@nodes, @element), $geometry)
        )
      )
    )
  ),
  array_find(@mean_distances, array_min(@mean_distances)) + 1
)

Github repo

Here, you’ll find a complete reproducible project with layers and expressions

Tutorial

One way to master QGIS field calculator expressions is to conceive them progressively, iteratively. Let’s decompose it.

Add a virtual field named nearest_with_nodes as integer which will store the index of the perfect target candidate, that is the index of the closest street line.

Open the Field Calculator window and follow this tutorial.

One street can be retrieved with the expression

get_feature_by_id('streets', 1)

Here it is the first object

To get the distance to this street from the current object, we calculate

distance(get_feature_by_id('streets', 1), $geometry)

But this distance is not ideal, because it is the closest space from one line to the other. If nodes are spread out, this does not reflect the « true » distance.

The « true » distance can be approached by calculating the distance from each of target object nodes to the current source line.

To get the nodes from the target line, we use

nodes_to_points(geometry(@parent) -- @nodes variable

To get the first node of the target line, we use

point_n(@nodes, 1)

To get the distance to this node

distance(point_n(@nodes, 1), $geometry)

But we have to compute it to each node.

We get the number of point with num_points(@nodes)

So we use generate_series to create the indexes of the nodes, going from 1 to the number of nodes. If 3 nodes, we’ll have [1,2,3]

generate_series(1, num_points(@nodes))

And we use array_foreach to calculate the distance to each node

array_foreach(
generate_series(1, num_points(@nodes)),
distance(point_n(@nodes, @element), $geometry)
)

Now we have an array of distances, let’s say [1, 3, 2], we want to get the average of the distance. For this, we use array_avg available after installing the ArrayPlus plugin.

array_avg(
  with_variable(
    'nodes',
    nodes_to_points(geometry(@parent)),
    array_foreach(
      generate_series(1, num_points(@nodes)),
      distance(point_n(@nodes, @element), $geometry)
    )
  )
)

And we calculate this mean distance from source line (pipeline) to each of the target lines (streets).

For this, we use aggregate function

aggregate(
  'streets',
  'array_agg',
  array_avg(
    with_variable(
      'nodes',
      nodes_to_points(geometry(@parent)),
      array_foreach(
        generate_series(1, num_points(@nodes)),
        distance(point_n(@nodes, @element), $geometry)
      )
    )
  )
)

Now, we have an array of mean distances, one mean distance per target line.

But we want to retrieve the min mean distance, and its index. Its index will be the one of the closest line.

To get the min distance, we use

array_min(@mean_distances)

We find the index of the min mean distance, in the array

array_find(@mean_distances, array_min(@mean_distances))

We add one because the first element of the array is indexed 0

array_find(@mean_distances, array_min(@mean_distances)) + 1

Here is the complete expression :

with_variable(
  'mean_distances',
  aggregate(
    'streets',
    'array_agg',
    array_avg(
      with_variable(
        'nodes',
        nodes_to_points(geometry(@parent)),
        array_foreach(
          generate_series(1, num_points(@nodes)),
          distance(point_n(@nodes, @element), $geometry)
        )
      )
    )
  ),
array_find(@mean_distances, array_min(@mean_distances)) + 1
)

So I hope you understood how to find the closest line to another one and the logic behind creating sophisticated expressions.

Display the street name

Let’s say the street layer has « name » attribute.

You can display as a label informations for the target line like this:

attribute(get_feature_by_id('streets', nearest_with_nodes), 'name')

Conclusion

Interestingly, you can see you can execute sophisticated tasks internally with expressions in the field calculator without making any SQLish virtual layers or applying magic python code.

That makes the Field Calculator a powerful and fun-to-use tool

To me, ArrayPlus should benefit from higher visibility as adding really interesting operators in QGIS. A user had to remain me of its existence. It would be nice if it joined the panoply of Field Calculator functions.

Bonus

What if you put many points on each target line ?

We can do this with interpolate_line_point

Here is the expression :

with_variable(
	'mult',
	10, 
	with_variable(
		'mean_distances',
		aggregate(
			'streets',
			'array_agg',
			with_variable(
				'feat',
				geometry(@parent),
				with_variable(
					'distances',
					array_foreach(
					generate_series(0, length(@feat), length(@feat)/((num_points(@feat)-1) * @mult)),
					distance($geometry, line_interpolate_point(@feat, @element))),
					array_avg(@distances)
				)
			)
		),
	array_find(@mean_distances, array_min(@mean_distances)) + 1
	)
)

geosparklines, a twitter map story

On the 10th of April, I created a map to see the evolution of COVID cases in France as spatialized sparklines.

I first saw it like a good use case on how to use some QGIS related functions (QGIS is an opensource geospatial tool)

« Les Artisans cartographes » mentioned the map on twitter. Thereafter, many people started retweeting it. I didn’t expect it to make so much « noise »

For me, this representation was an efficient way to see evolutions on a map because you could observe the dynamics at a single glance.

Maybe this instantaneous look can help in emergency cases because you don’t need to press play and wait for every frame to complete in order to get « the big picture ».

For me, there was nothing much innovative or disruptive about this map. It had been a long time the idea was in my head and I had already experimented it around census data.

I thought someone had already created something similar. But Alberto Cairo, a world-reknown dataviz expert, seemed to have never seen something like this before, which surprised me.

Kenneth Field, who is a reknown map maker compares this kind of map with Minard’s (a great honor !)

Read this wonderful article around Minard

Well before I put the QGIS code to create the map on github and explained it on this blog, some people on twitter started to create their own. Later, I saw some similar maps emerge.

That’s a gallery of these initiatives you’ll discover here.

Italy

Czech Republic

Also this one :

Holland

US

Wisconsin

Later, Kemper Smith made a US version

Washington Post

You can see that the Washington Post map is more sophisticated, with additional use of thickness to illustrate the current state (actual number of cases) and color to support the evolution

Note the Washington Post map was showcased on flowingdata, which is a reference media relaying dataviz experiments (also tutorials)

My experiments around this idea

An animated, worldwide and colored version

A hairy version

A scribbled version

A sparkline map of coviD-19 (or any name you’ll prefer)

This week-end, because of containment, I had lots of time to experiment around data and covid-19 cases, trying to figure out which graphical representation could draw the dynamics of the pandemic.

I saw a lot of graphics. The curves and their dynamics are crucial with the pandemic. Not only it is important to know what is the actual situation, but also how it was before, and how we can project ourselves to the future.

The above map is the map I made and shared on twitter :

It had quite a lot of success. It’s a mix of a line plot and a map. It gives a sense, at a single glance, of the dynamics, and looks good when animated.

After some feedback on twitter, I decied to improve the map with labels :

I only used opensource tools for it : qgis and imagemagick when animated. What’s great with opendata and open tools is that everyone can potentially integrate the data and experiment with no cost. You just have to let the ideas come, select them and find the ways to transcribe in the constrained environment of a machine : here, my computer. Somehow, the screen is my canvas painting, my mouse my brush and my keyboard my palette.

The goal of this article is to explain the method behind. It will give you a nice tutorial covering 3 qgis techniques :

  • 1-n relations
  • geometry generators
  • the use of arrays

Some next articles could get deeper into the subject with styling.

Download data

First and foremost, download the two datasets :

Here’s how the covid fr dataset looks like. It is in a long format with one row by day, unlike world John Hopskins worldwide dataset which is wide (with each column being a day)

Relation

We’ll add a relation between the departements dataset and the covid dataset. We use the department code INSEE_DEPT as a key connecting the parent table (departments layer) and the child table (covid dataset). Go to Project > Properties > Relation

Generative design

Next, go to the department layer and use a centroid symbol, but for this centroid symbol, you’ll use geometry generator to define the linestring.

As we want lines, choose line / polyline features :

And add this code in the expression (I’ll explain it just afterwards) :

with_variable(
'atlasid',
@atlas_featureid - 1,
with_variable(
'height',
100,
with_variable(
'width',
200000,
with_variable(
'nb_cas',
array_sort(
  relation_aggregate( 
    'donnees_ho_dep_DEPARTEMEN_INSEE_DEP', 
    'array_agg', 
     to_int("dc")+to_int("hosp")+to_int("rea")+to_int("rad"))
 ),
geom_from_wkt(
  'LINESTRING ('  || 
  x(centroid($geometry))||' '||y(centroid($geometry)) 
  ||','||
  array_to_string(
    array_foreach(
      generate_series(1, @atlasid), 
      -- X Coordinates
      to_string(x(centroid($geometry)) + 
      (@element/array_length(@nb_cas)*@width)) ||' '||
      -- Y coordinates
      to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
    ),
  ',') || ')'
)
))))

To understand this code, it’s important to split it into key parts. Obviously, I didn’t write it in a go. It was a trial-and-error process. Let me help you.

First, this line uses the relation we configured to store all the cases for each feature (here, each department), in an array

array_sort(
	relation_aggregate( 
		'donnees_ho_dep_DEPARTEMEN_INSEE_DEP', 
		'array_agg', 
		to_int("dc")+to_int("hosp")+to_int("rea")+to_int("rad")
	)
 )

We sort this array in ascending order of number of cases.

Note : you see the relation name being quite long, to select it, just go to the right panel and click on your relation name to put it in the code :

We store this array using the with_variable function in an array called nb_cas, to which we access with @nb_cas for future simplicity :

with_variable(
'nb_cas',
array_sort(
	relation_aggregate( 
		'donnees_ho_dep_DEPARTEMEN_INSEE_DEP', 
		'array_agg', 
		to_int("dc")+to_int("hosp")+to_int("rea")+to_int("rad")
	)
 )

Before constructing the linestring object, let’s consider this :

  • the height of the segment will be defined by the number of cases
  • and the spacing between nodes, simply, by the index of the node. One index = one node.
  • the line will start from the departement centroid

So, once we have this array of number of cases, we’ll iterate over it to get the consecutive nodes of the lines. We’ll just need the centroid coordinates, the indexes and the number of cases for each of them to build the lines

For the X coordinates of the nodes it’s this formula :

x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)

For the Y coordinates, it’s :

y(centroid($geometry)) + @nb_cas[@element]*@height

To get the index (@element), we use generate_series

generate_series(1, @atlasid), 
			-- X Coordinates
			to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
			-- Y coordinates
			to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
		)

To get the number of cases, we get the array element corresponding to the index :

@nb_cas[@element]

The iteration over indexes is made with array_for_each.

array_foreach(
			generate_series(1, @atlasid), 
			-- X Coordinates
			to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
			-- Y coordinates
			to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
		)

Once we have the array of coordinates, we paste them with array_to_string with the comma as delimiter, so that they can be interpreted as WKT string for a linestring

array_to_string(
		array_foreach(
			generate_series(1, @atlasid), 
			-- X Coordinates
			to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
			-- Y coordinates
			to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
		),
	',')

Our final geometry, transformed from a WKT string, will be :

geom_from_wkt(
	'LINESTRING ('  || 
	x(centroid($geometry))||' '||y(centroid($geometry)) ||','||
	array_to_string(
		array_foreach(
			generate_series(1, @atlasid), 
			-- X Coordinates
			to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
			-- Y coordinates
			to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
		),
	',') || ')'
)

As I said before, the first node here will be the centroid of the feature defined by its centroid. That’s why I pasted

x(centroid($geometry))||' '||y(centroid($geometry)) ||','||
	array_to_string(...

To make this the most programmatic as possible, I configured variables for the total width of the lines (@width), and the @height factor for each number of case. This way, it’s easy to change the appearance of the map

Also, I added a variable for the atlas feature id (@atlas_id) because I wanted to get the individual frames I could assemble as a GIF with ImageMagick.

For the animation, the trick consists in using an atlas and generating a series which will stop at the atlas current feature id, this one going from 0 to (n-1), n being the number of days.

So the first image will have one node (the centroid), the second one the centroid + the seconde node, the third image, the centroid + the two following nodes.

Note : why did I say 0 and (n-1) insteand of 1 and n ? It’s because the first element of an array, in qgis, is like in python : it has number 0.

That’s why my atlas_id variable has the value :

with_variable(
  'atlasid',
  @atlas_featureid - 1,
  @element
)

Maybe later I’ll write about styling the data underneath, but now, you have the mechanism and logic behind !

There’s a github repo 😉

You struggled to follow all these steps ? Wait ! I put all the project on a github repo with all the data you’ll need and a totally reproducible map. You’ll just have to follow the steps on the github page (actually, one single click on the project after downloadgin qgis)