Hands-on PostGIS, exploring the Geospatial capabilities
Tags: gis, postgres, node, geospatial
Published: 2020-07-03 10:54:22


I recently worked on a problem, the solution of which required geospatial querying and transformations. I initially thought of doing everything in-memory
using a JavaScript library called [turf](http://turfjs.org/). This was going fine until I used the actual dataset (I had around 24,000 Polygons having vertices ranging from 50 to 5,000), this is where I guess both the time and
space complexities of the operations I was doing via turf became too large, and it just gave up.

This is the point where I had to move to PostGIS in order to do those operations.

Someone might argue that out there is so much dedicated GIS software like QGIS, ArcGIS, why not use them. The answer is simple, it’s not only the GIS work that I am doing on my database, but there are also lots of regular CRUD operations as well, and PostGIS is just an extension to Postgres, which makes it perfect for this kind of use cases. And if you are thinking of suggesting me Mongo, I am gonna stop you right there, Mongo supports few geospatial queries but has no geo-transformation capabilities at all.

> In case you are worried to set up Postgres in replication mode and all, you can always use an [RDS](https://aws.amazon.com/rds/) instance. Recently I have been
hearing a lot of good things about [Amazon Aurora](https://aws.amazon.com/rds/aurora/) as well.

![](https://cdn-images-1.medium.com/max/1600/1*jSgZUQdYRBih_kbZAf9qKA.jpeg)
<span class="figcaption_hack">
But let’s define the problem statement first

### Problem Statement

Let’s say that you want to set up a fast-food chain in New York. But you want to
strategically place your restaurants to get the max customers, hence you decided
these criteria

1.  No restaurant should be near 205 meters of a KFC
1.  All restaurants should be in 300 meters of KFC
1.  All restaurants should be inside New York

You have to also generate the result in a GeoJSON file so that, that can be
drawn onto a map for better visualization.

> We are going to use GeoJSON to represent geometries throughout

Simple right? 😏 Let’s dive into code.

### Environment setup

1.  I am going to do this in [node.js](https://nodejs.org/en/), hence download and
install that.
1.  I am not going to install `postgres` and use a `docker image` instead, so
download and install [docker](https://www.docker.com/).
1.  Running these two commands in the image below consecutively will pull the docker
image for PostGIS and start a container.

![](https://cdn-images-1.medium.com/max/1600/1*CRshIjKRRMsu5Vv8diYtag.png)

> I am going to use `node.js` for this, don’t worry if you don’t understand
> `JavaScript/ TypeScript`, I am not going to use any language-specific `Object Relational Mapper` , I will execute raw PostGIS queries from code, which will make it language agnostic.

### Available Data

Let’s assume that we already have the boundary information of New York and all the KFCs there. For the sake of simplicity, I will not use the actual data but use a Polygon to represent NY and multiple Polygons to represent KFCs.

This is our imaginary NY (in gray), and KFCs (in red)

![](https://cdn-images-1.medium.com/max/1600/1*1kmlMzXIxxKH0xP_BXDOzQ.png)

These Polygon boundaries are available in `.json` files, which you can find
[here](https://github.com/dev-aritra/hands-on-postgis/tree/master/src/input).

### Project Setup

I am going to set up a very simple `node.js` + `TypeScript` console application.
If you are following along, you don’t have to again do that, you can download the template project from here. This already has all the `bootstrap` and `boilerplate` code ->[Branch_link](https://github.com/dev-aritra/hands-on-postgis/tree/step-1)

I am going to divide the entire problem into multiple small tasks and solve it step by step.

> I will share certain branches of my repo after each task, which contains the complete solution to that particular task. So don’t get baffled by small code snippet.
This is done so that you try to write the code yourself, which adds an active learning kind of flavor. Also, the complete solution would also be provided at the end of the article. Try to follow as much as you can, in case you fail, you
can always check out the branch for that task.

> If you didn’t understand this, then you probably need to learn **Git** first, which I will not explain, but you can learn it from [here](https://dev.to/unseenwizzard/learn-git-concepts-not-commands-4gjc).

### Task 1: Insert the NY Polygon and KFC Polygons into DB

Ok, so we need to first insert all the related data into the DB so that we can query/ operate on it.

To do so, I ended up writing these two files.

1.  The first is the `postgresPool.ts` file

```typescript
import { environment } from './../environment';
import * as pg from 'pg';

export async function getPostgresPool() {
    const pool = new pg.Pool({
        user: environment.postGisDBCreds.user,
        host: environment.postGisDBCreds.host,
        database: environment.postGisDBCreds.database,
        password: environment.postGisDBCreds.password,
        port: environment.postGisDBCreds.port
    });

    await pool.connect();
    return pool;
}
```

Which basically instantiates the `postgres connection pool`, which you can use to query the DB.

2. And the second is the `index.ts` file

```typescript
#!/usr/bin / env node

import * as pg from 'pg';
import * as newYork from './input/ny.json';
import * as kfcs from './input/kfcs.json';
import { getPostgresPool } from './dbClients/postgres-pool';
import { logGreen } from './logger';

main();

async function main() {
    const postgresPool = await getPostgresPool();
    await insertNYBoundary(postgresPool);
    await insertKFCBoundaries(postgresPool);
    logGreen('Execution complete');
}

async function insertNYBoundary(pgPool: pg.Pool) {
    const nyBoundaryTable = 'ny_boundary';
    const tableCreationQuery = buildTaleCreationQuery(nyBoundaryTable, pgPool);
    await pgPool.query(tableCreationQuery);

    const insertionQuery = buildInsertionQuery(JSON.stringify(newYork.geometry), nyBoundaryTable);
    await pgPool.query(insertionQuery);
}

async function insertKFCBoundaries(pgPool: pg.Pool) {
    const kfcBoundariesTable = 'kfc_boundaries';
    const tableCreationQuery = buildTaleCreationQuery(kfcBoundariesTable, pgPool);
    await pgPool.query(tableCreationQuery);

    const insertQueries = kfcs.features.map(kfc => buildInsertionQuery(JSON.stringify(kfc.geometry), kfcBoundariesTable));
    insertQueries.forEach(async insertQuery => await pgPool.query(insertQuery));
}

function buildTaleCreationQuery(tableName: String, pgPool: pg.Pool) {
    return 'create table if not exists ' + tableName + ' (id serial primary key, geom geometry)';
}

function buildInsertionQuery(data: any, table: String) {
    const queryPart1 = 'insert into ' + table + ' (geom) values (ST_GeomFromGeoJSON(\'';
    const queryPart3 = '\'));';
    const query = queryPart1.concat(data).concat(queryPart3);
    return query;
}
```


I know, I know it’s long, but let me explain, it’s pretty simple actually. The flow is like this

1.  Create a table with the name `ny_boundary`, having 2 columns, id, and geom.
1.  Read the geojson file `ny.json` and insert that into this table.
1.  Create a table with the name `kfc_boundaries`, having 2 columns, id, and geom.
1.  Read the geojson file `kfc.json` and insert that into this table.
1.  **buildTableCreationQuery** and **buildInsertionQuery** are basically 2 helper methods that generate the query for you given the data.

This is how the `create table` query would look if you run it as a SQL command

![](https://cdn-images-1.medium.com/max/1600/1*kIwAvSyX9UaCCJM3nJnTrQ.png)

And the `insert` SQL query

![](https://cdn-images-1.medium.com/max/1600/1*uCGv-BoOLF2OqDbQMNjgAQ.png)

We are using the *ST_GeomFromGeoJSON* function, cause the data we have as input is in geojson format.

And that concludes the task1, `woot woot` 🕺

And in case you could not follow along, here’s the ->
[branch_link](https://github.com/dev-aritra/hands-on-postgis/tree/step-1) that I promised. The boundary geojson files are at `/src/input`

### Task 2: Expand KFC boundaries by 205 meters and merge them if overlapping

Now this task has 2 sub-tasks.

1.  Expand the KFC boundaries by 205 meters, this will give us the area where we should not put our restaurants.
1.  We need to merge the expanded boundaries if there’s an overlap between any of them. Cause merged polygons look much better than overlapping polygons when rendered on a map. There are other benefits as well as it reduces the data size,
which can matter when dealing with huge data.

Ok, the code to achieve this
```typescript
async function expandBoundariesBy205Mtrs(pgPool: pg.Pool) {
    const level1Boundaries = 'level1_boundaries';
    const tableCreationQuery = buildTableCreationQuery(level1Boundaries, pgPool);
    await pgPool.query(tableCreationQuery);

    const expansionQuery = buildExpansionQuery(205);
    const expandedBoundaryInsertionQuery = 'insert into ' + level1Boundaries + ' (geom) ' + expansionQuery;
    await pgPool.query(expandedBoundaryInsertionQuery);
}

function buildExpansionQuery(distanceInMeters: number) {
    return 'select st_union(array (select st_buffer(geom::geography, ' + distanceInMeters + ' )::geometry from kfc_boundaries))'
}
```

Again, let me explain. I am doing two things.

1.  Creating a `level1_boundaries` table.
1.  Expanding all `kfc_boundaries` by 205 meters, merging them, and then inserting them into the `level1_boundaries` table.

I know the query part for the second operation might look a bit complex, a lot of things are going on there. So I will break down the parts of the query and try to explain what’s going on.

This is the query that we basically ran.
```sql
insert into level1_boundaries (geom) select st_union(array (select st_buffer(geom::geography, 205 )::geometry from kfc_boundaries));
```


     st_buffer(geom::geography, 205 )::geometry 
     kfc_boundaries

**st_buffer** -> This is the function which does the expansion operation, it takes two params, one is the geometry/ geography object, and the other is the radius for expansion. Which are `geom and 205` for our case.

**geom::geography **-> This is typecasting operation. We are taking the value of the geom column from kfc_boundaries and typecasting it to a `geography` object.
We need to this so so that the function considers the `geom` value as a EPSG:4326 geometry and hence it will consider the 205 as meters. If you don’t do this typecasting, the st_buffer function will consider the value 205 meters to be 205 degrees.

    ::geometry

We are again converting the result of the buffer operation to a `geometry` object, as `st_union` can only operate on `geometry` objects.

    st_union(array(...))

**st_union **-> This function merges the geometries returned from the st_buffer function if there are overlapping polygons it merges them into a single polygon, if the polygons are disjoint it creates a multipolygon out of them.

**array **-> As the kfc_boundaries table has multiple rows, the `select st_buffer ...` query will return an array, so to specify that it’s an array input, we are using this.

    insert into level1_boundaries (geom) ...

This basically inserts the result of the **st_union** operation into the `level1_boundaries` table.

Putting it together, this is how it looks now

![](https://cdn-images-1.medium.com/max/1600/1*QVVQLuxRkWqRDjoPqkLZWw.png)

*Blue polygons* -> KFC boundaries expanded by 205 meters

*Red Polygons* -> KFC boundaries

That’s the completion of task 2, and here’s the ->
[branch_link](https://github.com/dev-aritra/hands-on-postgis/tree/step-2)

### Task 3: Repeat step 2 but for a distance of 300 meters

Here’s the code

```typescript
async function expandBoundariesBy300Mtrs(pgPool: pg.Pool) {
    const level2Boundaries = 'level2_boundaries';
    const tableCreationQuery = buildTableCreationQuery(level2Boundaries, pgPool);
    await pgPool.query(tableCreationQuery);

    const expansionQuery = buildExpansionQuery(300);
    const expandedBoundaryInsertionQuery = 'insert into ' + level2Boundaries + ' (geom) ' + expansionQuery;
    await pgPool.query(expandedBoundaryInsertionQuery);
}
```

Nothing new, and here’s everything rendered on a map.

![](https://cdn-images-1.medium.com/max/1600/1*66IvjrYHyU_Lo_0x0vGP8A.png)

*Light green Polygons* -> KFC boundaries expanded by 300 meters

Here’s the ->
[branch_link](https://github.com/dev-aritra/hands-on-postgis/tree/step-3) for this task.

### Task 4: Subtract level 1 boundaries from level 2 boundaries to find out the
green zone

Now we have 2 MultiPolygons

1.  Boundary + 205 meters -> level 1
1.  Boundary + 300 meters -> level 2

We need to do

    level2 - level1

To find out the green zone, where we can set up our restaurants.

Here’s the code
```typescript
async function subtractLevel1FromLevel2(pgPool: pg.Pool) {
    const boundaryDifference = 'boundary_difference';
    const tableCreationQuery = buildTableCreationQuery(boundaryDifference, pgPool);
    await pgPool.query(tableCreationQuery);

    const level1 = (await pgPool.query('select geom from level1_boundaries')).rows[0].geom;
    const level2 = (await pgPool.query('select geom from level2_boundaries')).rows[0].geom;
    const query = "insert into " + boundaryDifference + " (geom) select ST_Difference(\'" + level2 + "\',\'" + level1 + "\');";
    await pgPool.query(query);
}
```


The flow is like this

1.  Create a table named `boundary_difference`
1.  Find the level1 boundary (the table has 1 row, hence for simplicity, I am using the first row only)
1.  Find the level2 boundary (the table has 1 row, hence for simplicity, I am using the first row only)
1.  Find the difference between these two objects using the `st_difference` function. It takes two geometries, finds out the difference between them, and
gives back the difference as a geometry object.
1.  Store the difference in `boundary_difference` table.

And here’s the outcome drawn on a map

![](https://cdn-images-1.medium.com/max/1600/1*80IOEFmD_tnXkuIA_H7NFg.png)

Here’s the ->
[branch_link](https://github.com/dev-aritra/hands-on-postgis/tree/step-4) for this task.

### Task 5: Find the intersection of NY boundary and green zones

You can see that some parts green zone which we just found out is going outside the NY boundary, and we do not want to set up our restaurants outside NY. Hence what we have to do now is to find out the intersection of the NY boundary and the green zone.

Here’s the code

```typescript
async function findIntersection(pgPool: pg.Pool) {
    const ny = (await pgPool.query('select geom from ny_boundary')).rows[0].geom;
    const difference = (await pgPool.query('select geom from boundary_difference')).rows[0].geom;
    const query = "select st_asgeojson(ST_Intersection(\'" + difference + "\',\'" + ny + "\'));";
    const newLocal = await pgPool.query(query);
    logCyan(JSON.stringify(newLocal.rows[0].st_asgeojson));
    writeFile('green_zones', JSON.stringify(newLocal.rows[0].st_asgeojson));
}
```

Flow

1.  Get the NY boundary and Difference calculated before
2.  Find out the intersection between them using `st_intersection`, which has a similar signature as `st_difference`, takes two geometries, and returns the
intersecting area as a geometry object.
3.  Convert the result into `geojson` using `st_asgeojson`
4.  Print the result in the console

And here’s the final picture, rendered on a map.

![](https://cdn-images-1.medium.com/max/1600/1*C6HLWI2p_mcZCDCENwAuqQ.png)

You can set-up you dream restaurant chain now 😆

And that my friend, concludes this tutorial on using some of the basic
functionalities of PostGIS to some pretty cool geospatial tasks.

You can download the complete solution from [here](https://github.com/dev-aritra/hands-on-postgis).

*****

Feel free to drop a comment if I got something wrong, I would love to improve.

If you learned anything, don’t forget to show support by pressing those like buttons.

You can find me on [Linkedin](https://www.linkedin.com/in/aritra-das-03843a8b/),
[GitHub](https://github.com/dev-aritra/). I would love to connect/ collaborate.

Happy learning! 😺