B.1.1 Attributes

The variable vector instantiated in the previous code lets us access the information stored in the shapefile. A general procedure is to go down the hierarchy of file > layer > feature to access different types of information. For a shapefile, one file only contains one layer, which in turn contains multiple features. Below are the lines of code for accessing the attributes associated with a feature.

We can write a short function (getattr.py) that returns a list of values in an attribute in a shapefile. We will use this function later in Section B.3.

B.1.2 Geometry and coordinates

We will now see how to retrieve the geometry and coordinates from the shapefile we just loaded. We only show how to get points of one polygon, but it is of course possible to get points for all polygons using a loop, as we will show later.

[Page 294]

B.1.3 Projecting points

We can transform a set of points in a geographic coordinate system (GCS) expressed in longitudes and latitudes into a projected coordinate system. Suppose we have a few locations from northern Cameroon that are recorded in WGS84 datum format, and we want to transform them into UTM zone 33N format. Here we will use the OSR (OGRSpatialReference) library, an OGR module for processing spatial reference systems, to transform between coordinate systems. OSR is based on the PROJ.4 Cartographic Projections Library.7 The general flow of work is to start by defining the source coordinate system, which is what we call gcs in the following code. Then we define the target coordinate system, called utmsr below. There are different ways to define a projected coordinate system, and we choose the easiest way using the EPSG code that is uniquely assigned to UTM 33N. With these two coordinate systems, we can define a transformation using the OSR function CoordinateTransformation, and here we name the transformation coordTransPcsUtm33. The known coordinates we have for the four locations in Cameroon are stored in a list of lists where the string in each element list refers to a place name. We first need to convert these coordinates into a point geometry in OGR and then we can simply use the Transform function to do the coordinate system transformation.

[Page 295]

B.1.4 Projecting and writing geospatial data

Now we turn to our shapefile that is in a GCS and this time we will project it to the Albers equal-area conic projection that is commonly used for the conterminous United States. Our shapefile comes with a .prj file and we can read the spatial reference information directly from there. The general idea here is to find each point and transform it to a new coordinate system and then write it to a ring, which is in turn added into a polygon feature. We save the transformed data into a new shapefile.

A general polygon in a shapefile may have holes, where each hole and the outline polygon are called rings. However, many counties in our data also have islands in a river or lake and these islands are not holes. In this case, they are organized as multipolygons. As the name suggests, each multipolygon includes a set of polygons, and each of those polygons may still contain holes. To streamline the process of handling such data, we design a specific function trans_rings that retrieves every ring in the input geometry that must be a polygon, which must not be a multipolygon but can have holes. (Having holes or not does not affect how this function works because holes and the outline polygon are all treated as rings and a polygon will have at least one ring.) This function returns a new polygon that contains transformed points. This polygon is then added into the feature for the new shapefile. For a feature that is only a simple polygon (with or without holes), we call this function once and get a new polygon. For a multipolygon feature, we will need to call this function multiple times, each time for a polygon in the multipolygon. Every time a new polygon is added, we create a new feature and add the polygon into the feature.

[Page 296]

[Page 297]

Figure B.1 Projecting spatial data. Left: original data in GCS. Right: data projected into Albers equal-area conic projection

The result is shown in Figure B.1. The two maps in this figure are drawn using Matplotlib using the following Python code (draw_shp_polygons.py).

[Page 298]

B.1.5 Adjacency matrix

Our final project using OGR is to create an adjacent matrix for a polygon shapefile. The procedure is simple: we continuously check each pair of polygons and test if they touch each other or one contains another. We use the OGR Touches and Contains functions to do so. To speed up the process, we design a function called env_touch to test if the bounding boxes (envelopes) of two polygons touch each other. We can write a function called geom_touch to test whether two OGR geometry objects touch each other by sharing a specified number of points (in parameter n0).

The adjacency_matrix function takes three input parameters: the path to the shapefile, the output format (M for matrix and L for list), and the number of shared points for polygons to be considered as adjacent (default is 1). In the test of the function, we use the US county shapefile and dump the final result in matrix form in a Python pickle file. In Python, pickle is a standard way of translating complicated objects so that we can store them in a text file and restore them for future use. The pickle file stored here is used in our main text (Section 9.2) for calculating Moran’s I.

[Page 299]

[Page 300]