Day 25

Thursday, May 4, 2017

Geographic Information Systems

Software specialized to deal with geospatial information.

Some applications

visualizing forest loss, other environmental changes
natural disaster relief management
tracking airline flights, package deliveries, your own hiking/jogging
population distributions and changes
election results, political divisions, recreational preferences
tectonic plate movement
weather and climate patterns
utility (water, sewer, gas, electricity, internet) network design, monitoring, and maintenance
self-driving cars
urban & regional planning, transportation
conducting warfare
adaptive farming

Shapefiles

Formerly proprietary format created by Esri, maker of ArcGIS.

A "shapefile" is actually a collection of files: .shp, .shx, .dbf, (.prj, .sbn, .sbx, ... ).

Let's take a quick look at the pieces and compare with Wikipedia description.

Coordinate systems and "datum"s

geodetic latitude depends on choice of datum (reference ellipsoid)

2017_05_04/20170504_105001.jpg

Before we use heavy-duty software ...

We can read shapefiles directly in Python using this shapefile module.

Exercise: draw census tracts

Unpack this shapefile of NY State US Census tracts, and use pylab.fill() to draw all the census tracts. Let pylab use its default coloring.

Things from shapefile module useful to know about:

foo = shapefile.Reader( filename )
foo.shapes()
shape.shapeType, shape.points, shape.parts
foo.records()
foo.fields

Let's use "Web Mercator" projection, i.e. Mercator projection

x = Rlon

y = Rlog(tan((π)/(4) + (lat)/(2))

except using geodetic latitude as if it were geocentric latitude.

2017_05_04/20170504_105012.jpg
import shapefile
for shape in f.shapes():
    #print(shape)
    #print(shape.shapeType)
    #print(shape.points)
    if len(shape.parts)>1: print(len(shape.points),shape.parts)
    #break
223 [0, 85]
1358 [0, 1213]
674 [0, 581]
805 [0, 412]
1276 [0, 945]
764 [0, 282]
416 [0, 334]
767 [0, 662]
1051 [0, 971]
760 [0, 396]
91 [0, 45]
560 [0, 434]
608 [0, 485, 493]
2340 [0, 1770, 1785, 1797, 1816, 1835, 1857, 1875, 1909, 1983, 2043, 2159, 2246]
362 [0, 245]
241 [0, 148, 214]
196 [0, 173]
315 [0, 149]
351 [0, 313]
333 [0, 169]
309 [0, 235]
564 [0, 299]
433 [0, 259]
210 [0, 122]
697 [0, 622]
334 [0, 279]
668 [0, 335]
1349 [0, 34, 94, 113, 187, 1233]
1457 [0, 12, 27, 1236, 1323, 1345, 1363]
878 [0, 19]
668 [0, 595]
929 [0, 588]
213 [0, 98]
319 [0, 46]
385 [0, 245]
1020 [0, 501]
284 [0, 257]
773 [0, 588]
1487 [0, 1262]
976 [0, 212]
123 [0, 8]
%pylab
R = 3959
f = shapefile.Reader('tl_2013_36_tract/tl_2013_36_tract.shp')
figure(figsize=(15,10))
for shape in f.shapes():
    #if len(shape.parts)>1: print(len(shape.points),shape.parts)
    #break
    a = array(shape.points)
    starts = shape.parts
    starts.append(len(shape.points))
    for i,part in enumerate(starts[:-1]):
        lon = a[starts[i]:starts[i+1],0]*(pi/180)
        lat = a[starts[i]:starts[i+1],1]*(pi/180)
        x = R*lon
        y = R*log(tan(pi/4+lat/2))
        fill(x,y,alpha=0.5)
tracts.png

Exercise: paint census tracts by data

Figure out how to get the dbf data into a pandas dataframe. Then redraw the census tracts using hexcolor() to color according to one or more of the data fields, such as "ALAND", "AWATER".

Note: to the casual observer, this looks similar to things we've done previously, but here we are drawing the picture from scratch (and we are reading shapefile formatted data).

QGIS - a spatial information viewer/editor

qgis_WienLogo.png

Let me know if you had any problems installing the Geographical Information System software QGIS.

Example data sources

Open Streetmap

openstreetmap.png

Location data from photos: get_lat_lon_exif_pil.py

Try this one:

20170503_184508.jpg

(I have location tagging turned on on my phone: take care if doing so yourself!!!)

location_tag.png

Location data from GPS phone app (MyTracks, GPSLogger): my long drive to work one day when I happened to be recording: Scajaquada_Expy.gpx

Location data from government sources: Abandoned buildings in Chicago from data.cityofchicago.org, real-time Chicago bus information that we've looked at previously.

Exercise: Use QGIS to plot your own point or track data

Use QGIS to take some list of lat-lon points generated or found by you and superimpose it on Open Streetmap. Upload a screenshot of your map to UBlearns.

Homework for next Thursday

Please read this paper, "You had me at hello: How phrasing affects memorability", which is the inspiration for the work that guest-presenter Prof. Rachael Hinkle will be showing us on Thursday, May 11.