Extracting MRT data with python

LTA Datamall provides MRT station data in .shp format. To extract this data in a python environment, I used pyshp. The .shp file contains the geo information, and the .dbf file contains more info about these points.

# import pyshp
import shapefile 

# read the file
myshp = open('data/TrainStations.shp', "rb")
mydbf = open('data/TrainStations.dbf', "rb")
sf = shapefile.Reader(shp=myshp, dbf=mydbf)
records = sf.shapeRecords()

# check what we have
print 'There are', len(records), 'shape objects in this file'

There are 101 shape objects in this file

Exploring the data

There are several possible types of shape. To find out what we have, we need to inspect the shape type. From pyshp's repo, we see the value 1 means these are of the type "POINT".

Using sf.fields(), we can see what fields are available in the .dbf file

print 'Type', sf.shapes()[0].shapeType
print sf.fields

Type 1
[('DeletionFlag', 'C', 1, 0), ['STN_NAM', 'C', 250, 0]]

Iterate over a few records to inspect the data

for record in records[:10]:
    print record.record[0], record.shape.points[0]

KRANJI [20083.391689386004, 45221.39571094048]
ONE NORTH [22889.0042042423, 31334.975294856424]
KALLANG [32228.671390099036, 32633.127592252626]
CHINESE GARDEN [16814.968460274642, 36033.31846462617]
TAMPINES [40446.578179615775, 37265.64583181421]
MACPHERSON [34337.17063360128, 34288.62835217322]
QUEENSTOWN [24964.715072883424, 30777.078770649703]
LITTLE INDIA [29779.6665958485, 32164.704880889523]
FARRER ROAD [25125.433166538987, 33302.965075892345]
EUNOS [35795.12505491873, 33560.10887875131]

Converting geodata to lat / lng

The coordinates look weird because Singapore has a special coordinate system called SVY21. Thankfully cgcai wrote an open source tool to convert these coordinates to lat/lng!

from utils.svy21 import SVY21
svy = SVY21()
for record in records[:10]:
    print svy.computeLatLon(record.shape.points[0][1], record.shape.points[0][0]) 
    

(1.4252389616569163, 103.76218030044807)
(1.2996558129131364, 103.78739365920227)
(1.3113959752999975, 103.87131493003949)
(1.3421443095143446, 103.73281406100547)
(1.353288524903338, 103.94515867441613)
(1.3263673703277745, 103.89026116371568)
(1.2946106656559997, 103.8060449033082)
(1.307159964327755, 103.84930939626489)
(1.317453904144435, 103.80748879117374)
(1.3197785646649076, 103.90336148486789)

Visualizing the data

Finally time to put this on a map!

import folium
sg_map = folium.Map(location=[1.38, 103.8], zoom_start=12)


for record in records:
    lat, lng = svy.computeLatLon(record.shape.points[0][1], record.shape.points[0][0])
    folium.RegularPolygonMarker(
        [lat, lng],
        popup=record.record[0].title(),
        fill_color='#0f0f0f',
        number_of_sides=4,
        radius=5
        ).add_to(sg_map)

sg_map