Spatial Data Processing with PySAL & Pandas

IPYNB

#by convention, we use these shorter two-letter names
import pysal as ps
import pandas as pd
import numpy as np

PySAL has two simple ways to read in data. But, first, you need to get the path from where your notebook is running on your computer to the place the data is. For example, to find where the notebook is running:

!pwd # on windows !cd

/Users/dani/code/gds_scipy16/content/part1

PySAL has a command that it uses to get the paths of its example datasets. Let's work with a commonly-used dataset first.

ps.examples.available()

['10740',
 'arcgis',
 'baltim',
 'book',
 'burkitt',
 'calemp',
 'chicago',
 'columbus',
 'desmith',
 'geodanet',
 'juvenile',
 'Line',
 'mexico',
 'nat',
 'networks',
 'newHaven',
 'Point',
 'Polygon',
 'sacramento2',
 'sids2',
 'snow_maps',
 'south',
 'stl',
 'street_net_pts',
 'taz',
 'us_income',
 'virginia',
 'wmat']

ps.examples.explain('us_income')

{'description': 'Per-capita income for the lower 47 US states 1929-2010',
 'explanation': [' * us48.shp: shapefile ',
  ' * us48.dbf: dbf for shapefile',
  ' * us48.shx: index for shapefile',
  ' * usjoin.csv: attribute data (comma delimited file)'],
 'name': 'us_income'}

csv_path = ps.examples.get_path('usjoin.csv')

f = ps.open(csv_path)
f.header[0:10]

['Name',
 'STATE_FIPS',
 '1929',
 '1930',
 '1931',
 '1932',
 '1933',
 '1934',
 '1935',
 '1936']

y2009 = f.by_col('2009')

y2009[0:10]

[32274, 32077, 31493, 40902, 40093, 52736, 40135, 36565, 33086, 30987]

Working with shapefiles

We can also work with local files outside the built-in examples.

To read in a shapefile, we will need the path to the file.

shp_path = '../data/texas.shp'
print(shp_path)

../data/texas.shp

Then, we open the file using the ps.open command:

f = ps.open(shp_path)

f is what we call a "file handle." That means that it only points to the data and provides ways to work with it. By itself, it does not read the whole dataset into memory. To see basic information about the file, we can use a few different methods.

For instance, the header of the file, which contains most of the metadata about the file:

f.header

{'BBOX Mmax': 0.0,
 'BBOX Mmin': 0.0,
 'BBOX Xmax': -93.50721740722656,
 'BBOX Xmin': -106.6495132446289,
 'BBOX Ymax': 36.49387741088867,
 'BBOX Ymin': 25.845197677612305,
 'BBOX Zmax': 0.0,
 'BBOX Zmin': 0.0,
 'File Code': 9994,
 'File Length': 49902,
 'Shape Type': 5,
 'Unused0': 0,
 'Unused1': 0,
 'Unused2': 0,
 'Unused3': 0,
 'Unused4': 0,
 'Version': 1000}

To actually read in the shapes from memory, you can use the following commands:

f.by_row(14) #gets the 14th shape from the file

<pysal.cg.shapes.Polygon at 0x10d8baa20>

all_polygons = f.read() #reads in all polygons from memory

len(all_polygons)

254

So, all 254 polygons have been read in from file. These are stored in PySAL shape objects, which can be used by PySAL and can be converted to other Python shape objects.

They typically have a few methods. So, since we've read in polygonal data, we can get some properties about the polygons. Let's just have a look at the first polygon:

all_polygons[0:5]

[<pysal.cg.shapes.Polygon at 0x10d8baba8>,
 <pysal.cg.shapes.Polygon at 0x10d8ba908>,
 <pysal.cg.shapes.Polygon at 0x10d8ba860>,
 <pysal.cg.shapes.Polygon at 0x10d8ba8d0>,
 <pysal.cg.shapes.Polygon at 0x10d8baa90>]

all_polygons[0].centroid #the centroid of the first polygon

(-100.27156110567945, 36.27508640938005)

all_polygons[0].area

0.23682222998468205

all_polygons[0].perimeter

1.9582821721538344

While in the Jupyter Notebook, you can examine what properties an object has by using the tab key.

polygon = all_polygons[0]

polygon. #press tab when the cursor is right after the dot

  File "<ipython-input-20-aa03438a2fa8>", line 1
    polygon. #press tab when the cursor is right after the dot
                                                              ^
SyntaxError: invalid syntax

Working with Data Tables

dbf_path = "../data/texas.dbf"
print(dbf_path)

../data/texas.dbf

When you're working with tables of data, like a csv or dbf, you can extract your data in the following way. Let's open the dbf file we got the path for above.

f = ps.open(dbf_path)

Just like with the shapefile, we can examine the header of the dbf file.

f.header

['NAME',
 'STATE_NAME',
 'STATE_FIPS',
 'CNTY_FIPS',
 'FIPS',
 'STFIPS',
 'COFIPS',
 'FIPSNO',
 'SOUTH',
 'HR60',
 'HR70',
 'HR80',
 'HR90',
 'HC60',
 'HC70',
 'HC80',
 'HC90',
 'PO60',
 'PO70',
 'PO80',
 'PO90',
 'RD60',
 'RD70',
 'RD80',
 'RD90',
 'PS60',
 'PS70',
 'PS80',
 'PS90',
 'UE60',
 'UE70',
 'UE80',
 'UE90',
 'DV60',
 'DV70',
 'DV80',
 'DV90',
 'MA60',
 'MA70',
 'MA80',
 'MA90',
 'POL60',
 'POL70',
 'POL80',
 'POL90',
 'DNL60',
 'DNL70',
 'DNL80',
 'DNL90',
 'MFIL59',
 'MFIL69',
 'MFIL79',
 'MFIL89',
 'FP59',
 'FP69',
 'FP79',
 'FP89',
 'BLK60',
 'BLK70',
 'BLK80',
 'BLK90',
 'GI59',
 'GI69',
 'GI79',
 'GI89',
 'FH60',
 'FH70',
 'FH80',
 'FH90']

So, the header is a list containing the names of all of the fields we can read. If we just wanted to grab the data of interest, HR90, we can use either by_col or by_col_array, depending on the format we want the resulting data in:

HR90 = f.by_col('HR90')
print(type(HR90).__name__, HR90[0:5])
HR90 = f.by_col_array('HR90')
print(type(HR90).__name__, HR90[0:5])

list [0.0, 0.0, 18.31166453, 0.0, 3.6517674554]
ndarray [[  0.        ]
 [  0.        ]
 [ 18.31166453]
 [  0.        ]
 [  3.65176746]]

As you can see, the by_col function returns a list of data, with no shape. It can only return one column at a time:

HRs = f.by_col('HR90', 'HR80')

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-25-1fef6a3c3a50> in <module>()
----> 1 HRs = f.by_col('HR90', 'HR80')


TypeError: __call__() takes 2 positional arguments but 3 were given

This error message is called a "traceback," as you see in the top right, and it usually provides feedback on why the previous command did not execute correctly. Here, you see that one-too-many arguments was provided to __call__, which tells us we cannot pass as many arguments as we did to by_col.

If you want to read in many columns at once and store them to an array, use by_col_array:

HRs = f.by_col_array('HR90', 'HR80')

HRs[0:10]

array([[  0.        ,   0.        ],
       [  0.        ,  10.50199538],
       [ 18.31166453,   5.10386362],
       [  0.        ,   0.        ],
       [  3.65176746,  10.4297038 ],
       [  0.        ,   0.        ],
       [  0.        ,  18.85369532],
       [  2.59514448,   6.33617194],
       [  0.        ,   0.        ],
       [  5.59753708,   6.0331825 ]])

It is best to use by_col_array on data of a single type. That is, if you read in a lot of columns, some of them numbers and some of them strings, all columns will get converted to the same datatype:

allcolumns = f.by_col_array(['NAME', 'STATE_NAME', 'HR90', 'HR80'])

allcolumns

array([['Lipscomb', 'Texas', '0.0', '0.0'],
       ['Sherman', 'Texas', '0.0', '10.501995379'],
       ['Dallam', 'Texas', '18.31166453', '5.1038636248'],
       ..., 
       ['Hidalgo', 'Texas', '7.3003167816', '8.2383277607'],
       ['Willacy', 'Texas', '5.6481219994', '7.6212251119'],
       ['Cameron', 'Texas', '12.302014455', '11.761321464']], 
      dtype='<U13')

Note that the numerical columns, HR90 & HR80 are now considered strings, since they show up with the single tickmarks around them, like '0.0'.

These methods work similarly for .csv files as well.

Using Pandas with PySAL

A new functionality added to PySAL recently allows you to work with shapefile/dbf pairs using Pandas. This optional extension is only turned on if you have Pandas installed. The extension is the ps.pdio module:

ps.pdio

<module 'pysal.contrib.pdutilities' from '/Users/dani/anaconda/envs/gds-scipy16/lib/python3.5/site-packages/pysal/contrib/pdutilities/__init__.py'>

To use it, you can read in shapefile/dbf pairs using the ps.pdio.read_files command.

shp_path = ps.examples.get_path('NAT.shp')
data_table = ps.pdio.read_files(shp_path)

This reads in the entire database table and adds a column to the end, called geometry, that stores the geometries read in from the shapefile.

Now, you can work with it like a standard pandas dataframe.

data_table.head()

The read_files function only works on shapefile/dbf pairs. If you need to read in data using CSVs, use pandas directly:

usjoin = pd.read_csv(csv_path)
#usjoin = ps.pdio.read_files(csv_path) #will not work, not a shp/dbf pair

usjoin.head()

The nice thing about working with pandas dataframes is that they have very powerful baked-in support for relational-style queries. By this, I mean that it is very easy to find things like:

The number of counties in each state:

data_table.groupby("STATE_NAME").size()

STATE_NAME
Alabama                  67
Arizona                  14
Arkansas                 75
California               58
Colorado                 63
Connecticut               8
Delaware                  3
District of Columbia      1
Florida                  67
Georgia                 159
Idaho                    44
Illinois                102
Indiana                  92
Iowa                     99
Kansas                  105
Kentucky                120
Louisiana                64
Maine                    16
Maryland                 24
Massachusetts            12
Michigan                 83
Minnesota                87
Mississippi              82
Missouri                115
Montana                  55
Nebraska                 93
Nevada                   17
New Hampshire            10
New Jersey               21
New Mexico               32
New York                 58
North Carolina          100
North Dakota             53
Ohio                     88
Oklahoma                 77
Oregon                   36
Pennsylvania             67
Rhode Island              5
South Carolina           46
South Dakota             66
Tennessee                95
Texas                   254
Utah                     29
Vermont                  14
Virginia                123
Washington               38
West Virginia            55
Wisconsin                70
Wyoming                  23
dtype: int64

Or, to get the rows of the table that are in Arizona, we can use the query function of the dataframe:

data_table.query('STATE_NAME == "Arizona"')

Behind the scenes, this uses a fast vectorized library, numexpr, to essentially do the following.

First, compare each row's STATE_NAME column to 'Arizona' and return True if the row matches:

data_table.STATE_NAME == 'Arizona'

0       False
1       False
2       False
3       False
4       False
5       False
6       False
7       False
8       False
9       False
10      False
11      False
12      False
13      False
14      False
15      False
16      False
17      False
18      False
19      False
20      False
21      False
22      False
23      False
24      False
25      False
26      False
27      False
28      False
29      False
        ...  
3055    False
3056    False
3057    False
3058    False
3059    False
3060    False
3061    False
3062    False
3063    False
3064    False
3065    False
3066    False
3067    False
3068    False
3069    False
3070    False
3071    False
3072    False
3073    False
3074    False
3075    False
3076    False
3077    False
3078    False
3079    False
3080     True
3081    False
3082    False
3083    False
3084    False
Name: STATE_NAME, dtype: bool

Then, use that to filter out rows where the condition is true:

data_table[data_table.STATE_NAME == 'Arizona']

We might need this behind the scenes knowledge when we want to chain together conditions, or when we need to do spatial queries.

This is because spatial queries are somewhat more complex. Let's say, for example, we want all of the counties in the US to the West of -121 longitude. We need a way to express that question. Ideally, we want something like:

SELECT
        *
FROM
        data_table
WHERE
        x_centroid < -121

So, let's refer to an arbitrary polygon in the the dataframe's geometry column as poly. The centroid of a PySAL polygon is stored as an (X,Y) pair, so the longitude is the first element of the pair, poly.centroid[0].

Then, applying this condition to each geometry, we get the same kind of filter we used above to grab only counties in Arizona:

data_table.geometry.apply(lambda x: x.centroid[0] < -121)\
                   .head()

0    False
1    False
2    False
3    False
4    False
Name: geometry, dtype: bool

If we use this as a filter on the table, we can get only the rows that match that condition, just like we did for the STATE_NAME query:

data_table[data_table.geometry.apply(lambda x: x.centroid[0] < -119)].head()

len(data_table[data_table.geometry.apply(lambda x: x.centroid[0] < -119)]) #how many west of -119?

109

Other types of spatial queries

Everybody knows the following statements are true:

1. If you head directly west from Reno, Nevada, you will shortly enter California.
2. San Diego is in California.

But what does this tell us about the location of San Diego relative to Reno?

Or for that matter, how many counties in California are to the east of Reno?

geom = data_table.query('(NAME == "Washoe") & (STATE_NAME == "Nevada")').geometry

lon,lat = geom.values[0].centroid

lon

-119.6555030699793

cal_counties = data_table.query('(STATE_NAME=="California")')

cal_counties[cal_counties.geometry.apply(lambda x: x.centroid[0] > lon)]

len(cal_counties)

58

This works on any type of spatial query.

For instance, if we wanted to find all of the counties that are within a threshold distance from an observation's centroid, we can do it in the following way.

But first, we need to handle distance calculations on the earth's surface.

from math import radians, sin, cos, sqrt, asin

def gcd(loc1, loc2, R=3961):
    """Great circle distance via Haversine formula

    Parameters
    ----------

    loc1: tuple (long, lat in decimal degrees)

    loc2: tuple (long, lat in decimal degrees)

    R: Radius of the earth (3961 miles, 6367 km)

    Returns
    -------
    great circle distance between loc1 and loc2 in units of R


    Notes
    ------
    Does not take into account non-spheroidal shape of the Earth



    >>> san_diego = -117.1611, 32.7157
    >>> austin = -97.7431, 30.2672
    >>> gcd(san_diego, austin)
    1155.474644164695


    """
    lon1, lat1 = loc1
    lon2, lat2 = loc2
    dLat = radians(lat2 - lat1)
    dLon = radians(lon2 - lon1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)

    a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
    c = 2*asin(sqrt(a))

    return R * c

def gcdm(loc1, loc2):
    return gcd(loc1, loc2)

def gcdk(loc1, loc2):
    return gcd(loc1, loc2, 6367 )

san_diego = -117.1611, 32.7157
austin = -97.7431, 30.2672
gcd(san_diego, austin)

1155.474644164695

gcdk(san_diego, austin)

1857.3357887898544

loc1 = (-117.1611, 0.0)
loc2 = (-118.1611, 0.0)
gcd(loc1, loc2)

69.13249167149539

loc1 = (-117.1611, 45.0)
loc2 = (-118.1611, 45.0)
gcd(loc1, loc2)

48.88374342930467

loc1 = (-117.1611, 89.0)
loc2 = (-118.1611, 89.0)
gcd(loc1, loc2)

1.2065130336642724

lats = range(0, 91)
onedeglon = [ gcd((-117.1611,lat),(-118.1611,lat)) for lat in lats]

import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(lats, onedeglon)
plt.ylabel('miles')
plt.xlabel('degree of latitude')
plt.title('Length of a degree of longitude')

<matplotlib.text.Text at 0x114174470>

san_diego = -117.1611, 32.7157
austin = -97.7431, 30.2672
gcd(san_diego, austin)

1155.474644164695

Now we can use our distance function to pose distance-related queries on our data table.

# Find all the counties with centroids within 50 miles of Austin
def near_target_point(polygon, target=austin, threshold=50):
    return gcd(polygon.centroid, target) < threshold 

data_table[data_table.geometry.apply(near_target_point)]

Moving in and out of the dataframe

Most things in PySAL will be explicit about what type their input should be. Most of the time, PySAL functions require either lists or arrays. This is why the file-handler methods are the default IO method in PySAL: the rest of the computational tools are built around their datatypes.

However, it is very easy to get the correct datatype from Pandas using the values and tolist commands.

tolist() will convert its entries to a list. But, it can only be called on individual columns (called Series in pandas documentation).

So, to turn the NAME column into a list:

data_table.NAME.tolist()[0:10]

['Lake of the Woods',
 'Ferry',
 'Stevens',
 'Okanogan',
 'Pend Oreille',
 'Boundary',
 'Lincoln',
 'Flathead',
 'Glacier',
 'Toole']

To extract many columns, you must select the columns you want and call their .values attribute.

If we were interested in grabbing all of the HR variables in the dataframe, we could first select those column names:

HRs = [col for col in data_table.columns if col.startswith('HR')]
HRs

['HR60', 'HR70', 'HR80', 'HR90']

We can use this to focus only on the columns we want:

data_table[HRs].head()

With this, calling .values gives an array containing all of the entries in this subset of the table:

data_table[['HR90', 'HR80']].values

array([[  0.        ,   8.85582713],
       [ 15.88562351,  17.20874204],
       [  6.46245315,   3.4507747 ],
       ..., 
       [  4.36732988,   5.2803488 ],
       [  3.72771194,   3.00003   ],
       [  2.04885495,   1.19474313]])

Using the PySAL pdio tools means that if you're comfortable with working in Pandas, you can continue to do so.

If you're more comfortable using Numpy or raw Python to do your data processing, PySAL's IO tools naturally support this.

Exercises

1. Find the county with the western most centroid that is within 1000 miles of Austin.
2. Find the distance between Austin and that centroid.