Click here to Skip to main content
13,771,011 members
Click here to Skip to main content
Add your own
alternative version

Tagged as

Stats

14K views
31 downloads
9 bookmarked
Posted 5 Sep 2014
Licenced CPOL

Creating New Map Regions in GeoJSON

, 5 Sep 2014
Rate this:
Please Sign up or sign in to vote.
Converting a PostGIS table to GeoJSON

Introduction

I sometimes need to show regions on a map. However, often there are no freely available (KML, shape, GeoJSON) geometries or topologies available that I can use, so I have to create them (almost) from scratch, e.g. based on existing city geometries. As these geometries tend to be too detailed, they have to be simplified as well.

As I tend to forget the way I need to do this, I thought it would be useful if I put it in a small article, so others (and myself) can use it as well. The example I'm using is based on Dutch open source data, so although the technique is generally applicable, the final result is probably not.

Background

In The Netherlands, the Dutch CENSUS (CBS) provides us with the geometries of all cities. New regions are often formed by creating a union of the existing city boundaries, and here I will demonstrate how to do that. You need to have a basic understanding of PostgreSQL and its spatial extension, PostGIS.

Method

Let me first show what the end result will look like for the 43 WMO regions (image created by GeoJSONLint).

WMO regions (end result)

Importing the CBS Data

  1. Download the data and extract the zip file.
  2. Using PgAdmin, create a new database (cbs) for the data.
  3. Enable the PostGis extension (in the new database, right-click on extensions and select 'New extension'. In the popup that opens, select postgis).
  4. Open a command prompt and load the data using the following command, where SRID = 28992 and DATABASE = cbs.
shp2pgsql -I -s <SRID> <PATH/TO/SHAPEFILE> | psql -d <DATABASE>

The final result is that we have the CBS shape files loaded into our database. In this case, we are only going to need the created gem_2013_v1 table. But before we can go on, we need to do one more thing: as the CBS data is in RD coordinates ('Rijksdriehoek', or SRID 28992), we need to set the spatial reference ID (SRID) for our geometries manually. I have no idea why shp2pgqsl does this automatically (am I missing something), but in order to support future SRID conversions, this is a requirement. So select the cbs table in PgAdmin, open a SQL window, and enter.

SELECT UpdateGeometrySRID('gem_2013_v1', 'geom', 28992);

BTW, you can check the SRID using the following command:

SELECT Find_SRID('cbs', 'gem_2013_v1', 'geom');

Creating a Mapping Between City Boundaries and the New Regions

As an example, I am going to create boundaries for the WMO health regions in The Netherlands. Although I can download it as a PDF or image, they unfortunately don't provide an option to download it in a GIS format, so I download the CSV format.

I'm not going to include an Excel tutorial, but I convert the downloaded data so I get something that looks like the following picture. Here, I create a match between a city code (gm_code) and a wmo_code.

region table

Next, I create a new table 'regios', right click on it and select import. Options: specify file, set to csv, specify it has a header, and the delimiter = ';'

Creating a New View with the Combined Geometry

Now it is time to create our new WMO geometries. Basically, this is just a 'simple' SQL query where you union the city boundaries, but as I am not very good in it, it took me some time to get it right. First, use a LEFT JOIN where you use the previously created regions table to match cities with wmo regions.

CREATE OR REPLACE VIEW public.wmo AS 
 SELECT r.wmo_naam,
    st_union(st_transform(g.geom, 4326)) AS geom,
    count(g.gm_code) AS aant_gem,
    sum(g.aant_inw) AS aant_inw,
    sum(g.aant_man) AS aant_man,
    sum(g.aant_vrouw) AS aant_vrouw,
    round(100::numeric * sum(g.p_00_14_jr * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_00_14_jr,
    round(100::numeric * sum(g.p_15_24_jr * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_15_24_jr,
    round(100::numeric * sum(g.p_25_44_jr * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_25_44_jr,
    round(100::numeric * sum(g.p_45_64_jr * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_45_64_jr,
    round(100::numeric * sum(g.p_65_eo_jr * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_65_eo_jr,
    round(100::numeric * sum(g.p_ongehuwd * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_ongehuwd,
    round(100::numeric * sum(g.p_gehuwd * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_gehuwd,
    round(100::numeric * sum(g.p_gescheid * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_gescheid,
    round(100::numeric * sum(g.p_verweduw * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_verweduw,
    round(100::numeric * sum(g.p_west_al * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_west_al,
    round(100::numeric * sum(g.p_n_w_al * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_n_w_al,
    round(100::numeric * sum(g.p_marokko * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_marokko,
    round(100::numeric * sum(g.p_ant_aru * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_ant_aru,
    round(100::numeric * sum(g.p_surinam * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_surinam,
    round(100::numeric * sum(g.p_turkije * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_turkije,
    round(100::numeric * sum(g.p_over_nw * g.aant_inw / 100::numeric) / sum(g.aant_inw)) AS p_over_nw,
    sum(g.auto_tot) AS auto_tot,
    sum(g.auto_hh) AS auto_hh,
    sum(g.auto_land) AS auto_land,
    sum(g.bedr_auto) AS bedr_auto,
    sum(g.motor_2w) AS motor_2w,
    sum(g.opp_tot) AS opp_tot,
    sum(g.opp_land) AS opp_land,
    sum(g.opp_water) AS opp_water
   FROM regios r
   LEFT JOIN gem_2013_v1 g ON r.gm_code::text = g.gm_code::text
  WHERE g.water::text = 'NEE'::text
  GROUP BY r.wmo_naam;

Specifically, note the following line, in which the geometries are unioned, and transformed to SRID 4326 (Web Mercator projection used by most web mapping engines):

st_union(st_transform(g.geom, 4326)) AS geom

Most other lines are just computing an average of existing properties.

Converting the Results to GeoJSON

GeoJSON is increasingly becoming more popular to display geometry data on the map, but it is only partially supported by PostGIS (i.e., you need to do some extra work to get it right). Fortunately, someone else already wrote a tutorial on this, and I only needed to change three things (marked in bold):

  • I limit the number of digits to 6
  • I've added the properties that I wished to export
  • I've specified the view that I was using

Run the following query in your SQL window to convert the output to GeoJSON.

SELECT row_to_json(fc)
 FROM (SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
 FROM (SELECT 'Feature' As type
    , ST_AsGeoJSON(lg.geom, 6)::json As geometry -- limit points to 6 digits
    , row_to_json((SELECT l FROM (SELECT aant_gem, aant_inw, aant_man, aant_vrouw, _
    p_00_14_jr, p_15_24_jr, p_25_44_jr, p_45_64_jr, p_65_eo_jr, p_ongehuwd, p_gehuwd, p_gescheid, _
    p_verweduw, p_west_al, p_n_w_al, p_marokko, p_ant_aru, p_surinam, p_turkije, p_over_nw, _
    auto_tot, auto_hh, auto_land, bedr_auto, motor_2w, opp_tot, opp_land, opp_water) As l
      )) As properties
   FROM wmo As lg   ) As f )  As fc;

As the resulting GeoJSON is very big, it cannot be displayed in the regular SQL window output, so I've selected to output the results to file using the menu 'Query|Export to file', where I specified the file name, no quoting and no column names. All being well, you should end up with a very big (15MB) GeoJSON file.

Simplifying the GeoJSON Data

Clearly, a 15MB GeoJSON file is not going to work on a SmartPhone, and also regular browsers will not be very thrilled about it. Fortunately, there is a tool that can help us simplify the GeoJSON (or Shape or TopoJSON) file, mapshaper. Although you can use it online for free, I was not very keen on uploading 15MB, so I've installed the local version using the Node Package manager.

npm install –g mapshaper 

And run it with:

?mapshaper-gui

It opens up a window allowing you to simplify your GeoJSON. Below, you see the before and after effect, and I could simplify the GeoJSON up to 0.34%.

MapShaper before and after simplifying

The finally generated GeoJSON can be verified using another free online tool, GeoJSONLint, which was used to create the image at the top.

The final GeoJSON file size was reduced from 15MB to a mere 95KB, so that worked very well.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Erik Vullings
Architect
Netherlands Netherlands
No Biography provided

You may also be interested in...

Pro

Comments and Discussions

 
QuestionYou can use ST_Simplify and ST_SimplifyPreserveTopology Pin
Omar.Pessoa12-Jan-18 2:07
memberOmar.Pessoa12-Jan-18 2:07 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Cookies | Terms of Use | Mobile
Web04-2016 | 2.8.181119.1 | Last Updated 5 Sep 2014
Article Copyright 2014 by Erik Vullings
Everything else Copyright © CodeProject, 1999-2018
Layout: fixed | fluid