TDM 40200 Project 13 - Introduction to Geospatial Data Analysis

While the map using raw count is useful for total volume, there can be bias incorporated; states such as Texas and California inherently will have more crashes as they are a lot more populated with larger land area and road and infrastructure counts. So to make comparison, we will normalize and create our visualization using density.

We will use 'ALAND' (land area) column we obtained from joining. It will first be converted to square km just for managing the numbers and readability.

Now we calculate the density measuring how concentrated the fatal crashes are in each state.

Please write the code to sort by either crash_cnt or crash_density column and make the highest values be shown first. Print the first 10 rows after.

The second print output shows:

We had Texas as the highest raw crash count previously, but here it disappeared from the top 10 list. Texas is geographically very big which makes crashes spread out over a large area. On the other hand, smaller states such as New Jersey and onwards were never seen and now they took over the top 10 list; they are geographically small but have a much denser infrastructure and high traffic volumes, so while the total crashes is lower, they can have higher concentration of incidents per square km. Lastly, District of Columbia appears as a bit of an outlier in density results due to its very small land area compared to other states; since the calculation used division of number of crashes by land area, even small number of crashes is amplified. This could lead to some sensitivity.

Try creating a mapped visualization for this, similar to what we did in question 3.

You will see that the map is generally same coloured with not much variation at all. This is not necessarily a data problem, but a scaling problem in visualization. Remember there was one data point significantly higher than the rest (DC); since a linear scale was used, all other values get compressed into the same bin for colouring. We can use norm=plt.matplotlib.colors.LogNorm(), which transforms into a log scale, to fix the visualization.

Now there are distinct shades across states that were not visible previously.

In previous questions, we looked at crash patterns and information at the state level, but there exists localized variations within a state. For instance, crashes in urban areas could differ from those occurring on rural places. To address this, we will examine pattern in a single state.

Local Focus

Let’s filter our data to Indiana only and look at the crash distribution. We will colour each crash point by whether it happened on a rural or urban road using the RUR_URB column (this indicates whether crash happened in a rural or an urban setting).

Below code will first tell us how many crashes happened in Indiana, and how they are split between rural and urban locations.

Now we visualize both on the same map.

We can also add a legend manually; to do this, we use:

Patches are geometric shapes we can draw on plot in matplotlib. This allows us to create patches and put in legend. In below code, we use in 'legend_elements' to define each parts and corresponding colours. We could create legends automatically, but this is another way to do it while having more control over placement and design, especially if we want to be specific in a more complex plot.

The code chunk above creates a list of legend elements. By using mpatches.patch(), Matplotlib knows to create coloured box in the legend. You can see the parameters where we can choose the legend colour box and label; we can also do custom formatting (it allows fstrings, value modifications, and getting certain stats like len()).

Some things we can notice include the way they cluster. Urban crashes seem to be more merged around major population and transportation areas, while rural crashes are more spread out across the state; also, a lot of them are most likely occurring near highways or longer distance roads with attributes such as higher speed limits.

Now let’s compare the fatality rates, and find out which of the two are more fatal on average. Please print the average fatality per crash for both areas.

You will see below statistic:

Rural crashes have higher average fatalities per incident compared to that of urban crashes. Taking real world scenario into account, we could infer they could be caused by factors such as higher speed limits as mentioned previously, less access to safety features and emergency responders, and varying road conditions.

Dissolve

Please print all unique values in the 'REGION' column of states_gdf.

US Census Bureau has geographic region definition that uses single integer value to represent four US regions: Northeast, Midwest, South, and West.

We will combine the crash stats with state geometry/boundaries. Below merge matches rows using 'NAME' (state name) and keeps all states in contiguous_states and add in the crash data.

Since coded values were used originally for each regions, we create readable region names. Below is a dictionary mapping numeric region codes to labels, where 1 = Northeast, 2 = Midwest, 3 = South, and 4 = West.

You will notice the the mapping operation is reflected. Take a look at REGION and REGION_NAME columns. For example, West Virginia has region 3 and the assigned name is South. You might also be able to infer some other information regarding how high or low crash counts and density are between regions relatively.

Now, we will use GeoPanda’s .dissolve() . A way to think of it is like performing groupby with shapes. As the name suggests, it merges geometries within a group into one feature and aggregates all rows in a group.

Print regions_gdf, and you will get:

Now, we are at the last step of normalizing and creating visualization; we will be able to compare regions visually. Please complete the code below. The process will be quite similar to what we did previously.

Then print the regional choropleth map, which will result in:

So far, we have analyzed the number of crashes that occurred in various areas; it may be a natural question to ask how far crashes are from a given location or object.

CRS for measurement

We have been using EPSG:4326 for coordinates in degrees, meaning we cannot compute accurate distances. For measurements, we will reproject the CRS to use meters. Here, we picked EPSG: 32616 (WGS 84 / UTM zone 16N), which is a valid choice for US midwestern states (contains Indiana and closer neighbouring states, and more specifically, approx. longitude 90°W and 84°W).

You will see the following output:

Clearly, the values changed from smaller real values representing decimal degrees, to larger values representing meters. These are distances in meters from the UTM zone; any calculations we do with these will be meter based results.

TIGER - PLACES

In the TIGER database, the cities and their boundaries are found under "Places", whose shapefiles contains data on cities, towns, etc, as well as Census Designated Places.

The Federal Information Processing Standard (FIPS) code for Indiana is 18, hence you see the dataset specificed with the code '18'. The code uniquely identifies state and its county level data.

Please print the shape, CRS, list of columns, and the head of the places dataset.

The CRS output will show 'EPSG:4269'. As we have already seen, this uses degrees and we cannot calculate the distance to a nearest point using this, so we will reproject to 'EPSG:32616' later on.

As well, the shape will be 'Shape: (696, 17)'; we have hundreds of places in the shapefile. Since this is too many for all to be used as reference points, we will filter for major cities.

Now, please write the code to pick top 10 cities largest by land area and save into selected_cities. Then, print only the 'NAME' and 'ALAND' of the selected cities.

Computing centroids for cities

Indianapolis, for example, is an irregular shape for city boundary. If we wanted to calculate distances, we would need a consistent, representative point. We can do so by using centroid, which is the geometric center of a polygon. This can be obtained automatically by GeoPandas. Before we do that, as previously mentioned, please convert the CRS to 'epsg=32616'.

Now, we compute centroids for the cities. Specifically, we can do so by using .geometry to access the geometry column, and .centroid to obtain the points representing the center of the region.

Please print the name and centroid in selected_cities.

Output:

Similarly as before, we no longer have latitude and longitude but values in meters. The first numbers is the distance from the central meridian of UTM zone, and the second number is the distance from the Equator.

Let’s create a visualization to see the centroids for each cities!

See below for more explanation of this section of the code:

We obtain the below map:

The blue polygons are TIGER provided city boundaries, and the red dots are computed centroids.

Distance - Crash & Nearest City

We will now explore distance calculation and sjoin_nearest() by using the centroids to see how far Indiana crashes occurred from nearest major city.

Print and inspect the head after this part.

Then, find which city appears most as the nearest city for crashes. In another word, how many times does each city appear in crashes_nearest?

We now create a spatial visualization to see how these distances are distributed across the state.

crashes_nearest had meters coordinates and not degrees, so we copy computed distances back onto the original dataframe.

In question 5, we plotted a fixed colour (one each for rural and urban), but here we will use ax.scatter(.. c=) to map distance which takes on continuous values.

The map will look like:

The map conveys multiple information simultaneously: geography, location of the crashes, and distances of crashes (another variable) from certain areas.