Skip to content

bag and postgis

Introduction

Bag delivers various products with horrible abbreviations, but luckily I found the translations here : visit

  • bag : basisregistratie adressen en gebouwen
  • bgt : basisregistratie grootschalige topografie
  • brt : basisregistratie topografie
  • dkk : digitale kadastrale kaart

Also it keeps mentioning TOP10, this seems to be under brt, which provides different types of postgis dumps :

More info on all the datasets is on the PDOK site : visit Nowhere does it mention what topXX means but my guess is TOPografie NL on on scale XX

By the way pdok means : Publieke Dienstverlening Op de Kaart (PDOK). Also look at their viewer : visit

  • top10NL : topografie Nederland schaal 1:10 ?
  • top50NL
  • top100NL
  • top250NL
  • top500NL
  • top1000NL : topografie NL schaal 1:1000 ?

This would mean top10 is the most detailed and thus largest dataset.

What we are after for klopt planner is BAG. To install there is a quick way and a slow way.

  • install a database dump from NLextract
  • install with nlextract software from the real sources at PDOK

We at least want to be able to perform the slow version to not be dependent on nlextract, but for fast updates, just try the pg_restore version.

The short way is not available for free anymore, use nlextract.

The next description is for when you buy a postgis database. We now extract it ourselves.

installation

The mongodb installation has been killed because of mongodb's policy change. The new implementation is for more than one reason postgis : visit

Since we don't have free downloads anymore, you will first have to generate the backup file as explained below (in nl extract), we now presume it is called bag-laatst.backup (since that was the old download name).

  • postgres was already a much tried and stable solution
  • postgres also has a GIS extension called postgis
  • the BAG actually comes as a postgis backup file at : visit

This method seems to fail after re-installation (maybe postgres/postgis version). Do complete install with nlextract (next chapter)

source : visit

  • You need to add a role called kademo
  • su - postgres
  • psql
  • create user kademo password 'kademo';
  • create database bagactueel;
  • \c bagactueel;
  • CREATE EXTENSION postgis;
  • pg_restore -v --dbname=bag bag-laatst.backup

Note that this creates a schema inside the bag database.

This means you cannot view the tables with dt :

changetitle
1
2
3
4
5
6
7
8
9
c bag;
d 
               List of relations
Schema |       Name        | Type  |  Owner   
--------+-------------------+-------+----------
public | geography_columns | view  | postgres
public | geometry_columns  | view  | postgres
public | spatial_ref_sys   | table | postgres
(3 rows)

These are just the tables created for the postgis extension. However the pg_restore command did install the 'bagactueel' schema :

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
SELECT schema_name FROM information_schema.schemata ;
schema_name     
--------------------
pg_toast
pg_temp_1
pg_toast_temp_1
pg_catalog
public
information_schema
bagactueel
(7 rows)

So you can view the tables need with commands like this :

changetitle
1
2
3
4
5
6
7
psql
c bag;
select * from bagactueel.pand limit 1;
# or 
dt bagactueel.*
# or even more ...
dt *.*

This also means that your source needs to connect to database 'bag' and do queries on tables : 'bagactueel.*'

If you want to play around with this use the bag-Amstelveen.backup file, which contains the bagaveen.* schema, so you can put both inside the bag database without clash.

Optionally you can list the tables use -l instead of -v Without the -v you get complete silence for several hours so at least to get some progress indication use -v.

If any errors get printed, drop the while database and try again.

  • drop database bag;
  • the list above.

If the pg_restore command keeps quiet you are on the right track. However this counts for postgres 9, on 12 this does not seems to lead to a usable database (no bagactueel table). So read next section to extract this manually into pg12.

nl extract

This program can read the BAG from source into various formats. The code is here:

visit

Manual in Dutch :

visit

This does seem to create a database without problems, but it is not usable. (no adres tabel for instance).

Here is the terse version, i created a working directory projects/klopt/bag. Make sure user postgres can write there. Also if you check that file out from bitbucket you don't need to install nlextract again.

The steps to do all this from scratch would have been :

changetitle
1
2
3
4
5
6
# as dev user (kees)
cd projects/klopt/bag
git clone https://github.com/nlextract/NLExtract.git
# now take some time making it python3 compatible !?!
# TODO : i never finished this and stepped back to python2 
# because the errors kept coming !1

Then to prepare the database, you need to install postgres+postgis and prepare pg_hba.conf to accept a connection.

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
vi /var/lib/postgres/data/pg_hba.conf 
# change all entries to md5 except postgres user itself
# and keep the ipv6 entry : host all all ::1/128 md5 !!
systemctl restart postgresql
su - postgres
createuser bag
createdb --owner bag -E UTF8 bag
psql bag
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
# change this next user/pass for production :
CREATE ROLE bag LOGIN ENCRYPTED PASSWORD 'bag';
CTRL-D # next command is shell again:

After this the bag database should contain 6 relations, 3 for each extension:

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
psql bag
d
             List of relations
  Schema  |       Name        |   Type   |  Owner   
----------+-------------------+----------+----------
public   | geography_columns | view     | postgres
public   | geometry_columns  | view     | postgres
public   | spatial_ref_sys   | table    | postgres
topology | layer             | table    | postgres
topology | topology          | table    | postgres
topology | topology_id_seq   | sequence | postgres
(6 rows)

The readme file says to create a separate postgis2 database on which you then install the extensions. After that you create the bag database with that as a template :

changetitle
1
2
# alternatively
createdb --owner bag -T postgis2 -E UTF8 bag

I have not tried that way, but next time see if it works, just always include the UTF8 encoding.

The database is prepared, now you need to prepare the extractor.

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
sudo su
apt-get install python-setuptools python-gdal unzip
# this one for when you get ERROR: fout nsmap in XML DOM processing
apt-get install python2-lxml
su - postgres 
cd ~kees/projects/klopt/bag
mkdir bag
cd bag
pip install --user psycopg2
pip install --user configparser
# get the latest data :
wget http://geodata.nationaalgeoregister.nl/inspireadressen/extract/inspireadressen.zip

Now edit the file NLExtract/bag/extract.conf, it is not an example, the script takes that file. Make sure you use the same username password in this file as in the db preparation. Then run this command, i run from the bag top directory so :

changetitle
1
2
3
python2 NLExtract/bag/src/bagextract.py -c
# or
NLExtract/bag/bin/bag-extract.sh -c

This prepares the database, to import you can use the -e option, which works on directories, zip files and xml files. So directory seems the easiest ?!:

changetitle
1
2
time NLExtract/bag/bin/bag-extract.sh -e inspireadressen.zip > log.txt
real    705m59.157s # 11 hours 45 minutes

Errors will scroll out of view when we don't redirect like this, of course if you want to see something happening you can tail -f the log file. Also the error are still printed on stderr so you can monitor what is going wrong easily. Note that this is by far the slower version with almost 12 hours.

To monitor the progress, here are the total number of xml files in each step :

  • NUM : 1149 about 60 minutes
  • LIG : 2
  • OPR : 34
  • PND : 1813 (3 sec * 1813 = 90 minutes)
  • STA : 5
  • VBO : 1636 (7 sec * 1636 = 190 min)
  • WPL : 1

total : real 422m23.663s is 7 hours

This is tested on the plain big zip file itself, but should work on the unzipped part files as well, and also on the .xml files that come out of that. This is certainly the easiest though but you may want to only do partly imports.

One issue may be the warning issued about self intersection:

changetitle
1
2
3
4
5
6
7
8
...
Warning 1: Self-intersection at or near point 45205.231345343811 397515.86791171448 0
Warning 1: Self-intersection at or near point 45205.231345343811 397515.86791171448 0
Warning 1: Self-intersection at or near point 38326.046951037657 399667.3494859084 0
Warning 1: Self-intersection at or near point 38326.046951037657 399667.3494859084 0
Warning 1: Self-intersection at or near point 50826.889600791386 399780.45922629244 0
Warning 1: Self-intersection at or near point 187614.36199999999 436143.283 0
...

See what problems these could cause .

Some more steps are needed for a complete database:

changetitle
1
2
3
4
5
6
python2 NLExtract/bag/src/bagextract.py -v -q NLExtract/bag/db/script/gemeente-provincie-tabel.sql 
time psql -d bag < NLExtract/bag/db/script/adres-tabel.sql
# 1507m41 (25 hours !!!)
cd NLExtract/bag/db/script/geocode
# this step still to be done (fails on psycopg2)
../../../bag/bin/bag-extract.sh -q geocode-tabellen.sql

Note that you will get errors on ST_Force_3D if you are using postgis3 because that function has been renamed ST_Force3D (no underscore !) So alter that if needed. It is in NLExtract/bag/db/script/adres-tabel.sql

After this step you will have these tables created :

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
d
                              List of relations
Schema |                     Name                     |   Type   |  Owner   
--------+----------------------------------------------+----------+----------
public | adres                                        | table    | bag
public | adres_gid_seq                                | sequence | bag
public | adresseerbaarobjectnevenadres                | table    | bag
public | adresseerbaarobjectnevenadres_gid_seq        | sequence | bag
public | adresseerbaarobjectnevenadresactueel         | view     | bag
public | adresseerbaarobjectnevenadresactueelbestaand | view     | bag
public | gemeente                                     | table    | bag
public | gemeente_gid_seq                             | sequence | bag
public | gemeente_woonplaats                          | table    | bag
public | gemeente_woonplaats_gid_seq                  | sequence | bag
public | gemeente_woonplaatsactueelbestaand           | view     | bag
public | geography_columns                            | view     | postgres
public | geometry_columns                             | view     | postgres
public | gt_pk_metadata                               | table    | bag
public | ligplaats                                    | table    | bag
public | ligplaats_gid_seq                            | sequence | bag
public | ligplaatsactueel                             | view     | bag
public | ligplaatsactueelbestaand                     | view     | bag
public | nlx_bag_info                                 | table    | bag
public | nlx_bag_info_gid_seq                         | sequence | bag
public | nlx_bag_log                                  | table    | bag
public | nlx_bag_log_gid_seq                          | sequence | bag
public | nummeraanduiding                             | table    | bag
public | nummeraanduiding_gid_seq                     | sequence | bag
public | nummeraanduidingactueel                      | view     | bag
public | nummeraanduidingactueelbestaand              | view     | bag
public | openbareruimte                               | table    | bag
public | openbareruimte_gid_seq                       | sequence | bag
public | openbareruimteactueel                        | view     | bag
public | openbareruimteactueelbestaand                | view     | bag
public | pand                                         | table    | bag
public | pand_gid_seq                                 | sequence | bag
public | pandactueel                                  | view     | bag
public | pandactueelbestaand                          | view     | bag
public | provincie                                    | table    | bag
public | provincie_gemeente                           | table    | bag
public | provincie_gemeente_gid_seq                   | sequence | bag
public | provincie_gemeenteactueelbestaand            | view     | bag
public | provincie_gid_seq                            | sequence | bag
public | raster_columns                               | view     | postgres
public | raster_overviews                             | view     | postgres
public | spatial_ref_sys                              | table    | postgres
public | standplaats                                  | table    | bag
public | standplaats_gid_seq                          | sequence | bag
public | standplaatsactueel                           | view     | bag
public | standplaatsactueelbestaand                   | view     | bag
public | verblijfsobject                              | table    | bag
public | verblijfsobject_gid_seq                      | sequence | bag
public | verblijfsobjectactueel                       | view     | bag
public | verblijfsobjectactueelbestaand               | view     | bag
public | verblijfsobjectgebruiksdoel                  | table    | bag
public | verblijfsobjectgebruiksdoel_gid_seq          | sequence | bag
public | verblijfsobjectgebruiksdoelactueel           | view     | bag
public | verblijfsobjectgebruiksdoelactueelbestaand   | view     | bag
public | verblijfsobjectpand                          | table    | bag
public | verblijfsobjectpand_gid_seq                  | sequence | bag
public | verblijfsobjectpandactueel                   | view     | bag
public | verblijfsobjectpandactueelbestaand           | view     | bag
public | woonplaats                                   | table    | bag
public | woonplaats_gid_seq                           | sequence | bag
public | woonplaatsactueel                            | view     | bag
public | woonplaatsactueelbestaand                    | view     | bag
(62 rows)

Still this needs some more indices for geocoding.

generation script

Will be described here ... This script will most probably take days to run as the test runs suggests.

Note that i want these conditions to be taken in to account:

  • all steps performed without human interference (download included)
  • error logging and progress logging to separate files
  • mail critical errors to kees@klopt.org
  • no interference with the current database, so import into a separate one or a separate system, better : incorporate date into the name, and switch in properties.conf
  • export the database by backup and restore it into the running systems.
  • also provide a manual work instruction.

Here is the history of commands i perform that should make up the script :

changetitle
1
2
3
4
5
# preliminary actions, only needed once 
git clone git@bitbucket.org:keesklopt/bag.git
cd bag
pip2 install -r requirements.txt
NLExtract/bag/bin/bag-extract.sh -c

Now every new update we probably want to keep running with the old database. The best scheme for that is to copy the database and fill the copy. But the best scenario is probably :

  • create a new db called bagimport
  • let all scripts run on that name, with extract.conf
  • when completed, rename bagimport to bagdate
  • switch the new name in /etc/properties.conf
  • restart klopt

That last step can also be automated, but it's probably best to just generate the new bag overnight/weekend and perform the switch and restart by hand.

So that means :

changetitle
1
2
3
4
5
6
sudo su 
cd /mnt/3/bag
wget http://geodata.nationaalgeoregister.nl/inspireadressen/extract/inspireadressen.zip
time NLExtract/bag/bin/bag-extract.sh -e inspireadressen.zip > log.txt
time python NLExtract/bag/src/bagextract.py -v -q NLExtract/bag/db/script/gemeente-provincie-tabel.sql
time psql -d bag -U bag < NLExtract/bag/db/script/adres-tabel.sql

The timed commands take (many) hours to complete. The script will be a weekend job !!

changetitle
1
NLExtract/bag/bin/bag-extract.sh -q NLExtract/bag/db/script/geocode/geocode-tabellen-full.sql

Note that the script say it is better to split the import data (bag_import) en geocode tables (bag_geocode) into schema's

A bit late.. so this first one is done in the public schema. In the script please do this division correctly !

But this will also get an error when you use postgres sql 11 or bigger. It will complain about lines like this :

changetitle
1
2
v.identificatie as adresseerbaarobject,
n.identificatie as nummeraanduiding,

The error is something like : column X is of type numeric but expression is of type character varying.

Altering all of these (3 lines) to read this will help :

CAST(v.identificatie as NUMERIC) as adresseerbaarobject,
CAST(n.identificatie as NUMERIC) as nummeraanduiding,

... # somewhat further down...
CAST(pv.identificatie as NUMERIC) as pandid,

Now this still gives lots of errors in the same area, but it does generate a usable database. So todo: fix these errors.

debian

The previous setup was on archlinux because the servert is archlinux, but for development, and also most production machines, it will be more likely to be debian machines. So here is the same install but now on debian buster.

I will however start with the already checked in git repo

changetitle
1
2
3
4
5
git clone git@bitbucket.org:keesklopt/bag.git
cd bag
git checkout devel # be sure to get the python2 version
wget http://geodata.nationaalgeoregister.nl/inspireadressen/extract/inspireadressen.zip
python -m pip install --user psycopg2 # when pip install does not work

This is a short version, but of course i already had plenty of tools installed so this did the trick...

the database

Note again that NLExtract did not deliver all tables as noted below, so we probably have to do a run with error checking, because for instance tables like openbareruimte are present but 'adres' is not. This would indicate a partial import.

The format of the database is like this

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
bag=# d
               List of relations
Schema |       Name        | Type  |  Owner   
--------+-------------------+-------+----------
public | geography_columns | view  | postgres
public | geometry_columns  | view  | postgres
public | raster_columns    | view  | postgres
public | raster_overviews  | view  | postgres
public | spatial_ref_sys   | table | postgres
(5 rows)

Note that you you d and NOT dt because dt only shows the tables.

If you view these than only geometry_columns and spatial_ref_sys are filled with data.

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
f_table_catalog | f_table_schema |          f_table_name          | f_geometry_column | coord_dimension | srid  |     type     
-----------------+----------------+--------------------------------+-------------------+-----------------+-------+--------------
bag             | bagactueel     | adres                          | geopunt           |               3 | 28992 | POINT
bag             | bagactueel     | adres_full                     | geopunt           |               3 | 28992 | POINT
bag             | bagactueel     | gemeente                       | geovlak           |               2 | 28992 | MULTIPOLYGON
BAG             | bagactueel     | ligplaats                      | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | ligplaatsactueel               | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | ligplaatsactueelbestaand       | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | pand                           | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | pandactueel                    | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | pandactueelbestaand            | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | provincie                      | geovlak           |               2 | 28992 | MULTIPOLYGON
bag             | bagactueel     | standplaats                    | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | standplaatsactueel             | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | standplaatsactueelbestaand     | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | verblijfsobject                | geopunt           |               3 | 28992 | POINT
bag             | bagactueel     | verblijfsobject                | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | verblijfsobjectactueel         | geopunt           |               3 | 28992 | POINT
bag             | bagactueel     | verblijfsobjectactueel         | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | verblijfsobjectactueelbestaand | geopunt           |               3 | 28992 | POINT
bag             | bagactueel     | verblijfsobjectactueelbestaand | geovlak           |               3 | 28992 | POLYGON
bag             | bagactueel     | woonplaats                     | geovlak           |               2 | 28992 | MULTIPOLYGON
bag             | bagactueel     | woonplaatsactueel              | geovlak           |               2 | 28992 | MULTIPOLYGON
bag             | bagactueel     | woonplaatsactueelbestaand      | geovlak           |               2 | 28992 | MULTIPOLYGON
(22 rows)

These should be familiar names if you worked with the java bag importer before. To see where the main data is you can use this postgres command :

changetitle
1
2
3
4
5
6
7
8
9
SELECT nspname || '.' || relname AS "relation",
 pg_size_pretty(pg_total_relation_size(C.oid)) AS "total_size"
FROM pg_class C
LEFT JOIN pg_namespace N ON (N.oid = C.relnamespace)
WHERE nspname NOT IN ('pg_catalog', 'information_schema')
   AND C.relkind <>'i'
   AND nspname !~ '^pg_toast'
ORDER BY pg_total_relation_size(C.oid) DESC
LIMIT 10;

Or just remove the limit at the end to see them all, the top list is :

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#                       relation                         | total_size 
---------------------------------------------------------+------------
 bagactueel.pand                                         | 8205 MB
 bagactueel.verblijfsobject                              | 5017 MB
 bagactueel.adres_full                                   | 4871 MB
 bagactueel.adres                                        | 3757 MB
 bagactueel.verblijfsobjectpand                          | 3568 MB
 bagactueel.verblijfsobjectgebruiksdoel                  | 2651 MB
 bagactueel.nummeraanduiding                             | 2507 MB
 bagactueel.openbareruimte                               | 75 MB
 bagactueel.woonplaats                                   | 31 MB
 bagactueel.adresseerbaarobjectnevenadres                | 20 MB
 bagactueel.standplaats                                  | 18 MB
 bagactueel.ligplaats                                    | 8432 kB
 bagactueel.gemeente                                     | 7976 kB
 public.spatial_ref_sys                                  | 4808 kB
 bagactueel.provincie                                    | 2896 kB
 bagactueel.nlx_bag_log                                  | 2296 kB
 bagactueel.gemeente_woonplaats                          | 1152 kB
 bagactueel.provincie_gemeente                           | 128 kB
 bagactueel.gt_pk_metadata                               | 24 kB
 bagactueel.nlx_bag_info                                 | 16 kB
 bagactueel.verblijfsobjectpand_gid_seq                  | 8192 bytes

This also show use different namespaces/schemas that are being used (or use dn;)

changetitle
1
2
3
4
5
6
7
dn;
List of schemas
   Name    |  Owner   
-----------+----------
bagactueel | postgres
public     | postgres
(2 rows)

Or if you want to print table data in all or a particular schema use :

changetitle
1
2
dt *.*;
dt bagactueel.*;

Clearly our data is in bagactueel. So all the queries will takes place there. These queries will also be used in the same format in the code so we can try to construct the ones we need on the command line.

find the coords of a postal code

In the BAG database postal code + housenumber are a complete identification of a location. So let's find my own coordinates : 3151aw 248. And then in reverse back to the address again + the written street and city data. I will traverse all tables to see where to get my data

adres

changetitle
1
select * from bagactueel.adres WHERE postcode='3151AW' AND huisnumm LIMIT 10;

Delivers :

changetitle
1
2
3
openbareruimtenaam  | huisnummer | huisletter | huisnummertoevoeging | postcode |  woonplaatsnaam  | gemeentenaam | provincienaam | typeadresseerbaarobject | adresseerbaarobject | nummeraanduiding | nevenadres |                              geopunt                               |                      textsearchable_adres                       |   gid   
---------------------+------------+------------+----------------------+----------+------------------+--------------+---------------+-------------------------+---------------------+------------------+------------+--------------------------------------------------------------------+-----------------------------------------------------------------+---------
 Prins Hendrikstraat |        248 |            |                      | 3151AW   | Hoek van Holland | Rotterdam    | Zuid-Holland  | VBO                     | 0599010000180840    | 0599200000414862 | f          | 01010000A04071000033333333AFC4F040C3F5285C1F181B410000000000000000 | '248':3 'hendrikstraat':2 'hoek':4 'holland':6 'prin':1 'van':5 | 5523581

So this already is a nice result, although the performance was less impressive. An index that i use now can be create with :

changetitle
1
CREATE INDEX adresses ON bagactueel.adres (postcode,huisnummer);

If it will be necessary to add annex(toevoeging) as well is not yet tested. It might slow down normal search !

This record contains everything we need :

  • streetname
  • house number
  • house letter
  • postcode
  • city
  • municipality
  • province
  • geopoint (coordinates)

The 'textsearchable_adres' is a column that is definitely worth a second look later on. !! The geopoint given is : 01010000A04071000033333333AFC4F040C3F5285C1F181B410000000000000000.

That's a bit unworkable, but if we can find the correct encoding it is resolvable into x and y coordinate.

For this you should read the explanations here : visit

If we use the standard function given there on this geometry :

changetitle
1
2
3
4
5
select ST_AsText(geopunt) from bagactueel.adres WHERE postcode='3151AW' AND huisnummer='248';
        st_astext            
--------------------------------
POINT Z (68682.95 443911.84 0)
(1 row)

This is clearly not decimal Lat Long but if we divide by 10000 we get coordinate 6.868295, 44.391185. Reverse these and paste into google maps :

changetitle
1
44.391185,6.868295

This puts us somewhere in France, so clearly we do not have the correct projection for this. The SRID for this is 28992 :

changetitle
1
SELECT srtext FROM spatial_ref_sys WHERE srid = 28992;

This results in pretty print :

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
PROJCS["Amersfoort / RD New",
  GEOGCS["Amersfoort",
     DATUM["Amersfoort",
        SPHEROID["Bessel 1841",6377397.155,299.1528128,
        AUTHORITY["EPSG","7004"]],
        TOWGS84[565.2369,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812],
        AUTHORITY["EPSG","6289"]],
     PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
     AUTHORITY["EPSG","4289"]],
  PROJECTION["Oblique_Stereographic"],
  PARAMETER["latitude_of_origin",52.15616055555555],
  PARAMETER["central_meridian",5.38763888888889],
  PARAMETER["scale_factor",0.9999079],
  PARAMETER["false_easting",155000],
  PARAMETER["false_northing",463000],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
  AXIS["X",EAST],
  AXIS["Y",NORTH],
  AUTHORITY["EPSG","28992"]]

See this :

visit

If you paste the original coordinates, 68682.95 and 443911.85 then hit transform you get these coordinates :

changetitle
1
51°58'36.838", 4°7'51.096"

Do not click show on the map, because while zooming in it loses position. Fill it in google maps search box, we have an exact match !!

This page has more information on the topic. It also seems we have to convert ourselves.

Also :

visit

I am ok with that because it's one linear step and probably the coordinate searching will be faster with the original encodings.

The coordinates we got were :

44.391185,6.868295

A number of solutions are presented to convert this :

  • qgis : apt-get install qgis
  • ogr2ogr : npm install ogr2ogr
  • proj4js : npm install proj4

The conversion we want is EPSG:28992 (Amersfoort RD/Nww) -> EPSG:4326 (WGS 84)

The parameters given earlier now can be used in libproj to get the right coordinates. A program found here provided the code, this is the original :

visit

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <iostream>
#include <proj_api.h>

using namespace std;

int main()
{
  double x = 138494.92605;
  double y = 467792.640021;

  char *srid28992 = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.04,49.91,465.84,-1.9848,1.7439,-9.0587,4.0772 +units=m +no_defs";
  char *srid4326 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";

  projPJ source = pj_init_plus(srid28992);
  projPJ target = pj_init_plus(srid4326);

  if(source==NULL || target==NULL)
    return false;

//  x *= DEG_TO_RAD;
//  y *= DEG_TO_RAD;

  int success = pj_transform(source, target, 1, 1, &x, &y, NULL );

  x *= RAD_TO_DEG;
  y *= RAD_TO_DEG;

  cout << success << endl << x << endl << y << endl;
}

However there are 2 things wrong with this example :

  • It has been forced into c++ for who knows why ?
  • The parameters for +towgs84 are wrong

First of all i removed all cout's and boolean's, then i altered the coordinates to the ones i am looking for. It then gave me a location somewhere close but a couple of streets away. Looking at the parameters i noticed they were slightly different from the output of postgis itself (see the pretty printed data earlier). Pasting those into place gave the correct result. The program now is :

changetitle
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include <stdio.h>
#include <proj_api.h>

int main()
{
  double x = 68682.95;
  double y = 443911.84;

  char *srid28992 = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs";
  char *srid4326 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";

  projPJ source = pj_init_plus(srid28992);
  projPJ target = pj_init_plus(srid4326);

  if(source==NULL || target==NULL)
    return 0;

//  x *= DEG_TO_RAD;
//  y *= DEG_TO_RAD;

  int success = pj_transform(source, target, 1, 1, &x, &y, NULL );

  x *= RAD_TO_DEG;
  y *= RAD_TO_DEG;

    printf("%f %fn", y, x);
}

The output :

changetitle
1
2
./a.out 
51.976899 4.130860

This program will be the basis for a conversion function in backend lib.

geocoding

Now the other way around, we not only want to search at this (exact) location but also have 'nearest point' functionality.

Since we have the coordinates in RDS format, we need the original coordinates for searching.

changetitle
1
2
double x = 68682.95;
double y = 443911.84;

A query that lists all addresses within a radius of 10 meters from the coordinates would be :

changetitle
1
2
3
4
5
6
7
SELECT postcode,huisnummer
FROM bagactueel.adres 
WHERE ST_DWithin(
     geom,
     ST_GeomFromText('POINT(68682.95 443911.84)',28992),
     10
   );

This renders 2 addresses, note that these are not necessarily in the order of closeness !

changetitle
1
2
3
4
5
postcode | huisnummer 
---------+------------
3151AW   |        246
3151AW   |        248
2 rows)

The performance of this operation is well within acceptable limits. It seems the necessary indices are already ok for this operation.

To get the distance as well and sort it use this query :

The performance of this operation is well within acceptable limits. It seems the necessary indices are already ok for this operation.