Plotting the Czech Cadastre Land Use with d3: Data Extraction (part I)

This post is the first part of the upcoming series summarizing the process of visualizing landuse data with bash, PostgreSQL and d3.js. Read other parts:

you’re reading it now
Plotting the Czech Cadastre Land Use with d3: Data Transformation (part II)
Plotting the Czech Cadastre Land Use with d3: Data Transformation (part III)

Czech Office for Surveying, Mapping and Cadastre has recently published lot of data via Atom feed. There’s pretty small and a bit boring dataset included, featuring quarterly updated landuse-related values for all 13,091 cadastral areas:

absolute number of land lots within given category (arable land, forests, etc.)
absolute area of land lots within given category

Data are published as CSV files linked from the Atom feed. Sadly, they come windows-1250 encoded, using Windows line endings, with trailing semicolons and header rows using diacritics.

ETL process

Before the d3 viz can be crafted, it’s necessary to:

extract CSV data from the URLs provided via the Atom feed
transform those data into a relational database, do some math
load data into a d3.js viz
profit (as usual)

Extract

#!/bin/bash
# extract.sh -f YYYYMMDD

while [[ $# -gt 1 ]]
do
key="$1"

case $key in
    -f|--file)
    FILE="$2"
    shift # past argument
    ;;
    *)
        # unknown option
    ;;
esac
shift # past argument or value
done

URL=http://services.cuzk.cz/sestavy/UHDP/UHDP-
CSVFILE=$FILE.csv
CSVUTF8FILE=${CSVFILE%.*}.utf.csv
URL+=$CSVFILE

echo "downloading $URL"
wget -q $URL -O $CSVFILE

if [[ $? != 0 ]]; then
    rm -f $CSVFILE
    echo "download failed"
    exit
fi

echo "converting to utf-8"
iconv -f WINDOWS-1250 -t UTF-8 $CSVFILE -o $CSVUTF8FILE &amp;&amp; \
echo "modifying ${FILE}"
sed -i 's/^M$//' $CSVUTF8FILE &amp;&amp; \
sed -i 's/\r$//' $CSVUTF8FILE &amp;&amp; \
sed -i 's/;*$//g' $CSVUTF8FILE &amp;&amp; \
sed -i '1d' $CSVUTF8FILE

echo "importing to database"
sed -e "s/\${DATE}/$FILE/g" extract.sql | psql -qAt --no-psqlrc

rm $CSVFILE $CSVUTF8FILE

This script downloads CSV file, deals with all the pitfalls mentioned above and, when done, copy command within extract.sql loads the data into a data_YYYYMMDD table. Putting all the files into the one table would have saved me a lot of transformation SQL, yet it didn’t feel quite right though.

Transform

See Plotting the Czech Cadastre Land Use with d3: Data Transformation (part II).

Load

See Plotting the Czech Cadastre Land Use with d3: Data Transformation (part III).

Introducing Blind Maps Project

I’d like to introduce you to my little pet project, which might just as well be awarded the first pet project I’ve ever completed, called Outline Maps of the World.

It’s a very simple, yet useful web application built on top of the great Leaflet library meant to help you get to know our world a bit better. As the name suggests, the app shows you, well… a blind map, and you try to fill as many features as you can.

The app is ready and can be used:

online at Outline Maps of the World with the map of your choice (if available)
offline, downloaded to your computer and filled with whatever data you want

What I find great about this project is the ease of adding new dataset. For starters, I filled it with data coming from Natural Earth:

CONUS states
European states
World capitals

If you wish, feel free to send me a pull request with GeoJSON data, I’ll be happy to have more datasets available! The process is described at the project homepage.

As you notice at the project homepage, there are two versions of the game available:

one lets you find map features by their names
the other one lets you type name highlighted feature (much tougher)

Have fun!

Degrees To Decimal With Javascript Reworked

Two years ago I was pretty happy with this little piece of code to transform degrees to the decimal value. Yesterday, I found a neater way to do the same:

let deg = [50, 30, 0];

function degToDec(prev, cur, curIndex) {
    return prev + cur / Math.pow(60, curIndex);
}

deg.reduce(degToDec);

Once you have an input array, that’s pretty much it. Love JavaScript.

Ogrinfo Output Formatting

Today my workmate asked if there was a way to see an attribute table other than importing spatial data into a PostGIS database. I told him about QGIS and while talking about other GIS stuff, I started thinking about pipes and how awesome it would be to actually format the output of the ogrinfo command.

Here it is. It is just a much longer way to do ogr2ogr -f "CSV" dest source, but sometimes you just have to experiment a bit.

#!/bin/bash
FILE=$1

function columns {
    ogrinfo $FILE -al -so | \
    sed '/Column/,$!d' | \
    sed '/Geometry Column/d' | \
    sed -e 's/Column =/\:/g' | \
    awk -F: '{print $1}' | \
    awk -v RS= -v OFS="|" '{$1 = $1} 1'
}

function data {
   ogrinfo $FILE -al | \
   sed '/OGRFeature/,$!d' | \
   sed '/POLYGON\|LINESTRING\|POINT/ d' | \
   sed -e 's/OGRFeature\(.*\)\://g' | \
   sed -e 's/.*\s*\(.*\)\s*=\s*//g' | \
   awk -v RS= -v OFS="|" '{$1 = $1} 1'
}

{ columns; data; }

The result can be piped to other bash functions, such as less or more. I call it ogrinfotable.

PostGIS Custom Function to Create Wind Rose

I’ve come across the beautiful GIS StackExchange question recently, asking how to draw a wind rose within PostGIS.

It’s pretty easy to accomplish this with a custom PLPGSQL procedure below, that takes line geometry, number of sections and radius of the inner circle as parameters.

CREATE OR REPLACE FUNCTION ST_WindRose(
    line geometry,
    directions int,
    radius numeric
)
RETURNS TABLE (
    id integer,
    geom geometry(LINESTRING)
)
AS $ST_WindRose$
BEGIN
    IF directions % 2 &lt;&gt; 0 THEN
        RAISE EXCEPTION 'Odd number of directions found, please provide even number of directions instead.';
    END IF;

IF radius &gt; ST_Length(line) THEN
    RAISE EXCEPTION 'Inner circle radius is bigger than the wind rose diameter, please make it smaller.';
END IF;

RETURN QUERY
WITH rose AS (
    SELECT
        ST_Rotate(_line, radians(360) / directions * dirs.id, ST_Centroid(_line)) _line
    FROM (
        SELECT line _line
    ) a
    CROSS JOIN (
        SELECT generate_series(1, directions / 2) id
    ) dirs
)
SELECT
    row_number() OVER ()::integer id,
    _line geom
FROM (
    SELECT _line FROM rose
    UNION ALL
    SELECT ST_ExteriorRing(ST_Buffer(ST_Centroid(line), radius, 30)) -- inner circle
    UNION ALL
    SELECT ST_ExteriorRing(ST_Buffer(ST_Centroid(line), ST_Length(line)/2, 30)) -- outer circle
) a;
END
$ST_WindRose$
LANGUAGE PLPGSQL;

Wind rose created with this function might look like the one below.

Run it as follows. The line parameter should be a simple straight line made of just two vertices.

SELECT * FROM ST_WindRose(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(0,1)), 12, 0.01);