Essays, Posts & Presentations

Population density map

Mads Bahrami

Author

Mads Bahrami

Title

Population density map

Description

Processing Facebook HDX data and make them computable in WL

Population Density Map

M. Bahrami

Wolfram Research, Inc.

Population density maps are among the most valuable datasets, e.g. when investigating the epidemic outbreak of an infectious disease. The Humanitarian Data Exchange (HDX) seems to have world’s most detailed population density datasets. These high-resolution maps, created by Facebook, estimate the number of people living within 30-meter grid tiles. In this short computational essay, I will discuss how to covert these maps into a computable datasets in the Wolfram Language, with desirable size and shape of grid tiles. An an example, I shall focus on Qatar and show how to get desired maps using GeoHistogram, a very powerful built-in function in the Wolfram Language.

The HDX high resolution population density maps can be obtained through their website. The dataset is provided in in two formats, GeoTIFF and CSV, for distribution of various populations. I will download and process the CSV file of the overall population density. The same methods can be applied for any other desired population dataset. Since Qatar’s dataset is small, I shall import it directly. But for larger datasets, I suggest downloading the file and then importing it.

Import zipped data from HDX and get the CSV file:

In[]:=

qatar=Import

,"population_qat_2019-07-01.csv";

Return the number of tiles on the original data:

In[]:=

Length[qatar]

Out[]=

307646

Return the headings of the original data:

In[]:=

qatar〚1〛

Out[]=

{Lat,Lon,Population}

Generate a list of two sublists as geo positions and a corresponding population:

In[]:=

{geo,ppl}={qatar〚2;;,;;2〛,qatar〚2;;,-1〛};

Generate a list whose elements are {geoposition->population}:

In[]:=

data=Thread[(GeoPosition/@geo)ppl];

Let’s look at the data using the default bins and counts of GeoHistogram.

Generate a density histogram:

In[]:=

GeoHistogram[data,PlotStyleEdgeForm[Gray],GeoScaleBar"Metric",PlotLegendsAutomatic,GeoRange

100

]

Out[]=

Using GeoHistogram, one can set the size and shape of the grid tiles.

Specify the diameter of each tile:

In[]:=

Table[GeoHistogram[data,i*

,PlotStyleEdgeForm[Gray],GeoScaleBar"Metric",PlotLegendsAutomatic,GeoRange

100

,PlotLabel"Grid Tile ("<>ToString[i]<>"km)"],{i,10,25,5}]

Out[]=

In[]:=

Table[GeoHistogram[data,i*

^2,PlotStyleEdgeForm[Gray],GeoScaleBar"Metric",PlotLegendsAutomatic,GeoRange

100

,PlotLabel"Grid Tile ("<>ToString[i]<>"

)"],{i,100,800,200}]

Out[]=

One can also set three different tile shapes: “Hexagon”, “Rectangle” and “Triangle”.

In[]:=

Table[GeoHistogram[data,{i,250*

^2},PlotStyleEdgeForm[Gray],GeoScaleBar"Metric",PlotLegendsAutomatic,GeoRange

100

,PlotLabel"Tile Shape: "<>i],{i,{"Hexagon","Rectangle","Triangle"}}]

Out[]=

Therefore, given a tile shape and size, we need to extract the geo positions of tiles generated by GeoHistogram and their corresponding values.

In[]:=

gh=GeoHistogram[data,{"Hexagon",250*

^2}];

If you look at the full form of the GeoHistogram output, you will notice that the grid tiles are wrapped in Tooltip.

In[]:=

Short[polygons=Cases[gh,Tooltip[Polygon[y_?MatrixQ],x_]Polygon[GeoPosition[y]]x,∞]]

Out[]//Short=

PolygonGeoPosition

Number of points: 6

Lat bounds: {51.4,51.6}

Lon bounds: {26.,26.2}

1.28428×

,PolygonGeoPosition

Number of points: 6

Lat bounds: {51.4,51.6}

Lon bounds: {26.2,26.4}

174644.,63,PolygonGeoPosition

Number of points: 6

Lat bounds: {50.9,51.1}

Lon bounds: {27.,27.1}

102.292

However, the actual positions of polygon points are not correct.

In[]:=

GeoGraphicsEdgeForm[Red],polygons[[All,1]],EdgeForm[Black],Polygon

Qatar

COUNTRY



Out[]=

A closer look shows that the polygon points in the GeoHistogram output are also projected using GeoGridPosition. Accordingly, in order to extract polygons with correct geo positions, we need to include the projection using GeoGridPosition too.

In[]:=

polygonCounts=Cases[gh,Tooltip[y:Polygon[_?MatrixQ],x_]:>GeoPosition[GeoGridPosition[y,(GeoProjection/.Options[gh,GeoProjection])]]x,∞];

In[]:=

GeoGraphicsEdgeForm[Red],polygonCounts[[All,1]],EdgeForm[Black],Polygon

Qatar

COUNTRY



Out[]=

Given this result, we now put the data in the format of a list with sublists as {polygon points, center of polygon, population density} and export it as CSV file (might be useful for users of other programming languages).

In[]:=

Export[NotebookDirectory[]<>"QatarPopulationDensity.csv",Transpose[{polygonCounts〚All,1,1,1〛,(RegionCentroid/@polygonCounts〚All,1〛)〚All,1〛,polygonCounts〚All,2〛}],TableHeadings{"Polygon points {Lat,Long}","Center of Polygon","Population Density"}]

Out[]=

/Users/mohammadb/Nextcloud/Mads/Exploration/COVID19/Epidemic/QatarPopulationDensity.csv

Summary

Given the data from HDX, generate grid tiles of population density for a given shape and size of tiles:

In[]:=

PopulationDensityGridTiles[rawData_,shape_/;MemberQ[{"Hexagon","Rectangle","Triangle"},shape],size_]:=Module[{gh,geo,ppl,dt},{geo,ppl}={rawData〚2;;,;;2〛,rawData〚2;;,-1〛};dt=Thread[(GeoPosition/@geo)ppl];gh=GeoHistogram[dt,{shape,size}];Cases[gh,Tooltip[y:Polygon[_?MatrixQ],x_]:>GeoPosition[GeoGridPosition[y,(GeoProjection/.Options[gh,GeoProjection])]]x,∞]]

Define colors, bins and legends for the GeoGraphics function:

In[]:=

colors=ResourceFunction["ColorBrewerData"]["Reds",9];bins={5,10,20,30,40,50,100,500}*1000;cfunc[x_]:=Which@@Riffle[Join[x≤#&/@bins,{x>Last[bins]}],colors]legend=SwatchLegend[colors,ToString/@Join[{"<"<>ToString@First[bins]},ToString@#1<>"-"<>ToString@#2&@@@Partition[bins,2,1],{">"<>ToString@Last[bins]}],LegendMarkersGraphics[{EdgeForm[Black],Opacity[1],Rectangle[]}],LegendLabel"Overall population",LegendFunction(Framed[#,RoundingRadius5]&),LegendMargins5];

Visualize each tile with a color specified by cfunc, defined before:

In[]:=

LegendedGeoGraphicsTooltip[{EdgeForm[Gray],GeoStyling[Opacity[.8],cfunc[#2]],#1},#2]&@@@PopulationDensityGridTiles[qatar,"Rectangle",{

}],EdgeForm[Black],FaceForm[],Polygon

Qatar

COUNTRY

,GeoScaleBar{"Metric","Imperial"},GeoRange

,legend

Out[]=

Acknowledgment

I’d like to thank my colleagues, Alan Joyce (Wolfram|Alpha team) and José Martín-García (Algorithms R&D team), for their valuable input in helping with datasets and important suggestions on code.

Cite this as: Mads Bahrami, "Population density map" from the Notebook Archive (2020), https://notebookarchive.org/2020-04-dushoma