Geospatial Analytics using Python and SQL

This adapted from Purdue’s HONR 39900-054: Foundations in Geospatial Analytics. For more information and examples, please head to the course’s GitHub site: github.com/gouldju1/honr39900-foundations-of-geospatial-analytics.

*In[5]:*

*In[21]:*

*In[6]:*

*Out[6]:*

*In[7]:*

*Out[7]:*

We have a LOT to cover today…here’s a snapshot:

Textbook Chapter/Section: 11.1.1
Textbook Start Page: 310

In many cases, you may have a city where records are broken out by districts, neighborhoods, boroughs, or precincts because you often need to view or report on each neighborhood separately. Sometimes, however, for reporting purposes you need to view the city as a single unit. In this case you can use the ST_UNION aggregate function to amass one single multipolygon from constituent multipolygons.

For example, the largest city in the United States is New York, made up of the five storied boroughs of Manhattan, Bronx, Queens, Brooklyn, and Staten Island. To aggregate New York, you first need to create a boroughs table with five records—one multipolygon for each of the boroughs with littorals:

Then you can use the ST_Union spatial aggregate function to group all the boroughs into a single city, as follows:

*In[10]:*

example_11_1_1 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="city") #Note the change in geom_col!

*In[11]:*

*Out[11]:*

*In[12]:*

Let’s work through an example in the San Francisco area using our table of cities. This example lists cities that straddle multiple records, how many polygons each city straddles, and how many polygons you’ll be left with after dissolving boundaries within each city.

*In[13]:*

ex_11_1_1_2 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="dummy") #Note the change in geom_col!

*In[14]:*

*Out[14]:*

From the code above, you know that ten cities have multiple records, but you’ll only be able to dissolve the boundaries of Brisbane and San Francisco, because only these two have fewer polygons per geometry than what you started out with.

Below, you will aggregate and insert the aggregated records into a table called ch11.distinct_cities. You then add a primary key to each city to ensure that you have exactly one record per city…

*In[15]:*

*In[16]:*

*In[18]:*

*In[19]:*

*In[20]:*

*Out[20]:*

*In[21]:*

In the above code, you create and populate a new table called ch11.distinct_cities. You use the ST_Multi function to ensure that all the resulting geometries will be multipolygons and not polygons. If a geometry has a single polygon, ST_Multi will upgrade it to be a multipolygon with a single polygon. Then you cast the geometry using typmod to ensure that the geometry type and spatial reference system are correctly registered in the geometry_columns view. For good measure we also put in a primary key and a spatial index.

Textbook Chapter/Section: 11.1.2
Textbook Start Page: 312

In the past two decades, the use of GPS devices has gone mainstream. GPS Samaritans spend their leisure time visiting points of interest (POI), taking GPS readings, and sharing their adventures via the web. Common venues have included local taverns, eateries, fishing holes, and filling stations with the lowest prices. A common follow-up task after gathering the raw positions of the POIs is to connect them to form an unbroken course.

In this exercise, you’ll use Australian track points to create linestrings. These track points consist of GPS readings taken during a span of about ten hours from afternoon to early morning on a wintry July day. We have no idea of what the readings represent. Let’s say a zoologist fastened a GPS around the neck of a roo and tracked her for an evening. The readings came in every ten seconds or so, but instead of creating one line-string with more than two thousand points, you’ll divide the readings into 15-minute intervals and create separate linestrings for each of the intervals.

ST_MakeLine is the spatial aggregate function that takes a set of points and forms a linestring out of them. You can add an ORDER BY clause to aggregate functions; this is particularly useful when you need to control the order in which aggregation occurs. In this example, you’ll order by the input time of the readings.

*In[22]:*

example_11_1_2 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="geom")

*In[23]:*

*Out[23]:*

*In[24]:*

First you create a column called track_period specifying quarter-hour slots starting on the hour, 15 minutes past, 30 minutes past, and 45 minutes past. You allocate each GPS point into the slots and create separate linestrings from each time slot via the GROUP BY clause. Not all time slots need to have points, and some slots may have a single point. If a slot is devoid of points, it won’t be part of the output. If a slot only has one point, it’s removed. For the allocation, you use the data_part function and the modulo operator.

Within each slot, you create a linestring using ST_MakeLine. You want the line to follow the timing of the measurements, so you add an ORDER BY clause to ST_MakeLine.

The SELECT inserts directly into a new table called aussie_run. (If you’re not running this code for the first time, you’ll need to drop the aussie_run table first.) Finally, you query aussie_run to find the number of points in each linestring using ST_NPoints, subtracting the time of the last point from the time of the first point to get a duration, and using ST_Length to compute the distance covered between the first and last points within the 15-minute slot. Note that you cast the geometry in longitude and latitude to geography to ensure you have a measurable unit—meters.

Textbook Chapter/Section: 11.2.1
Textbook Start Page: 314

Clipping uses one geometry to cut another during intersection.Clipping uses one geometry to cut another during intersection. In this section, we’ll explore other functions available to you for clipping and splitting.In this section, we’ll explore other functions available to you for clipping and splitting.

As the name implies, clipping is the act of removing unwanted sections of a geometry, leaving behind only what’s of interest. Think of clipping coupons from a newspaper, clipping hair from someone’s head, or the moon clipping the sun in a solar eclipse.

Difference and symmetric difference are operations closely related to intersection. They both serve to return the remainder of an intersection. ST_Difference is a non—commutative function, whereas ST_SymDifference is, as the name implies, commutative.As the name implies, clipping is the act of removing unwanted sections of a geometry, leaving behind only what’s of interest. Think of clipping coupons from a newspaper, clipping hair from someone’s head, or the moon clipping the sun in a solar eclipse.

Difference and symmetric difference are operations closely related to intersection. They both serve to return the remainder of an intersection. ST_Difference is a non—commutative function, whereas ST_SymDifference is, as the name implies, commutative.

Difference functions return the geometry of what’s left out when two geometries intersect. When given geometries A and B, ST_Difference(A,B) returns the portion of A that’s not shared with B, whereas ST_SymDifference(A,B) returns the portion of A and B that’s not shared. Difference functions return the geometry of what’s left out when two geometries intersect. When given geometries A and B, ST_Difference(A,B) returns the portion of A that’s not shared with B, whereas ST_SymDifference(A,B) returns the portion of A and B that’s not shared.

Textbook Chapter/Section: 11.2.2
Textbook Start Page: 316

You just learned that using a linestring to slice a polygon with ST_Difference doesn’t work. For that, PostGIS offers another function called ST_Split. The ST_Split function can only be used with single geometries, not collections, and the blade you use to cut has to be one dimension lower than what you’re cutting up.

The following listing demonstrates the use of ST_Split.

*In[35]:*

*Out[35]:*

*In[36]:*

The ST_Split(A,B) function always returns a geometry collection consisting of all parts of geometry A that result from splitting it with geometry B, even when the result is a single geometry.

Because of the inconvenience of geometry collections, you’ll often see ST_Split combined with ST_Dump, as in the above listing, or with ST_CollectionExtract to simplify down to a single geometry where possible.

Textbook Chapter/Section: 11.2.3
Textbook Start Page: 318

Dividing your polygon into regions using shapes such as rectangles, hexagons, and triangles is called tessellating. It’s often desirable to divide your regions into areas that have equal area or population for easier statistical analysis. In this section, we’ll demonstrate techniques for achieving equal-area regions.

Tessellating looks like this:

Textbook Chapter/Section: 11.2.3
Textbook Start Page: 327

One common reason for breaking a geometry into smaller bits is for performance. The reason is that operations such as intersections and intersects work much faster on geometries with fewer points or smaller area. If you have such a need, the fastest sharding function is the ST_Subdivide(postgis.net/docs/ST_Subdivide.html) function. ST_Subdivide is a set returning function that returns a set of geometries where no geometry has more than max_vertices designated. It can work with areal, linear, and point geometries.

This next example divides Queens such that no polygon has more than 1/4th total vertices of original.

*In[44]:*

example_11_2_3_3 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="geom") #Note change of geom_col! example_11_2_3_3.head()

*Out[44]:*

*In[45]:*

Here is a color-coded example of the above plot:

You will often find ST_Subdivide paired with ST_Segmentize. ST_Segmentize is used to even out the segments so that when ST_Subdivide does the work of segmentizing, the areas have a better chance of being equal.

In the above queries with ST_Subdivide use LATERAL join constructs. In the case of set-returning functions, the LATERAL keyword is optional so you will find people often ommitting LATERAL even though they are doing a LATERAL join.

Textbook Chapter/Section: 11.3.1
Textbook Start Page: 330

There are several reasons why you might want to break up a linestring into segments: - To improve the use of spatial indexes—a smaller linestring will have a smaller bounding box. - To prevent linestrings from stretching beyond one unit measure. - As a step toward topology and routing to determine shared edges.

If you have long linestrings where the vertices are fairly far apart, you can inject intermediary points using the ST_Segmentize function. ST_Segmentize adds points to your linestring to make sure that no individual segments in the linestring exceed a given length. ST_Segmentize exists for both geometry and geography types.

For the geometry version of ST_Segmentize, the measurement specified with the max length is in the units of the spatial reference system. For geography, the units are always meters.

In the following listing, you’re segmentizing a 4 vertex linestring into 10,000-meter segments.

*In[46]:*

example_11_3_1_1 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="geom") #Note change of geom_col! example_11_3_1_1.head()

*Out[46]:*

*In[47]:*

In this listing, you start with a 4-point linestring. After segmentizing, you end up with a 3,585-point linestring where the distance between any two adjacent points is no more than 10,000 meters. You cast the geography object to a geometry to use the ST_NPoints function. The ST_NPoints function does not exist for geography.

As mentioned earlier, ST_Subdivide is often used with ST_Segmentize. This is because ST_Subdivide cuts at vertices, so if you wanted to then break apart your linestring into separate linestrings, you can use ST_Subdivide, but need to first cast to geometry and then back to geography as follows:

*In[49]:*

example_11_3_1_2 = gpd.GeoDataFrame.from_postgis(sql=sql, con=engine, geom_col="geo") #Note change of geom_col! example_11_3_1_2.head()

*Out[49]:*

*In[50]:*

Textbook Chapter/Section: 11.4.2
Textbook Start Page: 337

The scaling family of functions comes in four overloads, one for 2D ST_Scale(geometry, xfactor, yfactor), one for 3D ST_Scale(geometry, xfactor, yfactor, zfactor), one for any geometry dimension but will scales the coordinates based on a factor specified as a point ST_Scale(geometry, factor), and the newest addition introduced in PostGIS 2.5 is a version that will scale same amount but about a specfied point ST_Scale(geom, factor, geometry origin). (postgis.net/docs/ST_Scale.html).

Scaling takes a coordinate and multiplies it by the factor parameters. If you pass in a factor between 1 and -1, you shrink the geometry. If you pass in negative factors, the geometry will flip in addition to scaling. Below shows an example of scaling a hexagon.

*In[68]:*

*Out[68]:*

*In[69]:*

*In[70]:*

*In[71]:*

You start with a hexagonal polygon #1 and shrink and expand the geometry in the X and Y directions from 50% of its size to twice its size by using a cross join that generates numbers from 0 to 2 in X and 0 to 2 in Y, incrementing .5 for each step #2. The results are shown below:

The scaling multiplies the coordinates. Because the hexagon starts at the origin, all scaled geometries still have their bases at the origin. Normally when you scale, you want to keep the centroid constant. See below for an exampe:

*In[72]:*

Here you scale a hexagon from half to twice its size in the X and Y directions about the centroid so the centroid remains unaltered. See below:

Textbook Chapter/Section: 11.4.3
Textbook Start Page: 340

ST_RotateX, ST_RotateY, ST_RotateZ, and ST_Rotate rotate a geometry about the X, Y, or Z axis in radian units. ST_Rotate and ST_RotateZ are the same because the default axis of rotation is Z. ST_RotateX, ST_RotateY, ST_RotateZ, and ST_Rotate rotate a geometry about the X, Y, or Z axis in radian units. ST_Rotate and ST_RotateZ are the same because the default axis of rotation is Z.

These functions are rarely used in isolation because their default behavior is to rotate the geometry about the (0,0) origin rather than about the centroid. You can pass in an optional point argument called pointOrigin. When this argument is specified, rotation is about that point; otherwise, rotation is about the origin.

The following listing rotates a hexagon about its centroid in increments of 45 degrees.

*In[78]:*

*In[80]:*

The above plot looks like a big blob…here is a better representation: