<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ibis and Geospatial Operations
One of the most popular extensions to PostgreSQL is PostGIS,
which adds support for storing geospatial geometries,
as well as functionality for reasoning about and performing operations on those geometries.
This is a demo showing how to assemble ibis expressions for a PostGIS-enabled database.
We will be using a database that has been loaded with an Open Street Map
extract for Southern California.
This extract can be found here,
and loaded into PostGIS using a tool like ogr2ogr.
Preparation
We first need to set up a demonstration database and load it with the sample data.
If you have Docker installed, you can download and launch a PostGIS database with the following
Step1:  Next, we download our OSM extract (about 400 MB)
Step2:  Finally, we load it into the database using ogr2ogr (this may take some time)
Step3:  Connecting to the database
We first make the relevant imports, and connect to the PostGIS database
Step4:  Let's look at the tables available in the database
Step5:  As you can see, this Open Street Map extract stores its data according to the geometry type.
Let's grab references to the polygon and line tables
Step6:  Querying the data
We query the polygons table for shapes with an administrative level of 8,
which corresponds to municipalities.
We also reproject some of the column names so we don't have a name collision later.
Step7:  We can assemble a specific query for the city of Los Angeles,
and execute it to get the geometry of the city.
This will be useful later when reasoning about other geospatial relationships in the LA area
Step8:  Let's also extract freeways from the lines table,
which are indicated by the value 'motorway' in the highway column
Step9:  Making a spatial join
Let's test a spatial join by selecting all the highways that intersect the city of Los Angeles,
or one if its neighbors.
We begin by assembling an expression for Los Angeles and its neighbors.
We consider a city to be a neighbor if it has any point of intersection
(by this critereon we also get Los Angeles itself).
We can pass in the city geometry that we selected above when making our query by marking it as a literal value in ibis
Step10:  Now we join the neighbors expression with the freeways expression,
on the condition that the highways intersect any of the city geometries
Step11:  Combining the results
Now that we have made a number of queries and joins, let's combine them into a single plot.
To make the plot a bit nicer, we can also load some shapefiles for the coast and land
<ASSISTANT_TASK:>
Python Code:
# Launch the postgis container.
# This may take a bit of time if it needs to download the image.
!docker run -d -p 5432:5432 --name postgis-db -e POSTGRES_PASSWORD=supersecret mdillon/postgis:9.6-alpine
Explanation: Ibis and Geospatial Operations
One of the most popular extensions to PostgreSQL is PostGIS,
which adds support for storing geospatial geometries,
as well as functionality for reasoning about and performing operations on those geometries.
This is a demo showing how to assemble ibis expressions for a PostGIS-enabled database.
We will be using a database that has been loaded with an Open Street Map
extract for Southern California.
This extract can be found here,
and loaded into PostGIS using a tool like ogr2ogr.
Preparation
We first need to set up a demonstration database and load it with the sample data.
If you have Docker installed, you can download and launch a PostGIS database with the following:
End of explanation
!wget https://download.geofabrik.de/north-america/us/california/socal-latest.osm.pbf
Explanation: Next, we download our OSM extract (about 400 MB):
End of explanation
!ogr2ogr -f PostgreSQL PG:"dbname='postgres' user='postgres' password='supersecret' port=5432 host='localhost'" -lco OVERWRITE=yes --config PG_USE_COPY YES socal-latest.osm.pbf
Explanation: Finally, we load it into the database using ogr2ogr (this may take some time):
End of explanation
import os
import geopandas
import ibis
%matplotlib inline
client = ibis.postgres.connect(
    url='postgres://postgres:supersecret@localhost:5432/postgres'
)
Explanation: Connecting to the database
We first make the relevant imports, and connect to the PostGIS database:
End of explanation
client.list_tables()
Explanation: Let's look at the tables available in the database:
End of explanation
polygons = client.table('multipolygons')
lines = client.table('lines')
Explanation: As you can see, this Open Street Map extract stores its data according to the geometry type.
Let's grab references to the polygon and line tables:
End of explanation
cities = polygons[polygons.admin_level == '8']
cities = cities[
    cities.name.name('city_name'),
    cities.wkb_geometry.name('city_geometry')
]
Explanation: Querying the data
We query the polygons table for shapes with an administrative level of 8,
which corresponds to municipalities.
We also reproject some of the column names so we don't have a name collision later.
End of explanation
los_angeles = cities[cities.city_name == 'Los Angeles']
la_city = los_angeles.execute()
la_city_geom = la_city.iloc[0].city_geometry
la_city_geom
Explanation: We can assemble a specific query for the city of Los Angeles,
and execute it to get the geometry of the city.
This will be useful later when reasoning about other geospatial relationships in the LA area:
End of explanation
highways = lines[(lines.highway == 'motorway')]
highways = highways[
    highways.name.name('highway_name'),
    highways.wkb_geometry.name('highway_geometry')
]
Explanation: Let's also extract freeways from the lines table,
which are indicated by the value 'motorway' in the highway column:
End of explanation
la_neighbors_expr = cities[
    cities.city_geometry.intersects(
        ibis.literal(la_city_geom, type='multipolygon;4326:geometry')
    )
]
la_neighbors = la_neighbors_expr.execute().dropna()
la_neighbors
Explanation: Making a spatial join
Let's test a spatial join by selecting all the highways that intersect the city of Los Angeles,
or one if its neighbors.
We begin by assembling an expression for Los Angeles and its neighbors.
We consider a city to be a neighbor if it has any point of intersection
(by this critereon we also get Los Angeles itself).
We can pass in the city geometry that we selected above when making our query by marking it as a literal value in ibis:
End of explanation
la_highways_expr = highways.inner_join(
    la_neighbors_expr,
    highways.highway_geometry.intersects(la_neighbors_expr.city_geometry)
).materialize()
la_highways = la_highways_expr.execute()
la_highways.plot()
Explanation: Now we join the neighbors expression with the freeways expression,
on the condition that the highways intersect any of the city geometries:
End of explanation
ocean = geopandas.read_file('https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_ocean.zip')
land = geopandas.read_file('https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip')
ax = la_neighbors.dropna().plot(figsize=(16,16), cmap='rainbow', alpha=0.9)
ax.set_autoscale_on(False)
ax.set_axis_off()
land.plot(ax=ax, color='tan', alpha=0.4)
ax = ocean.plot(ax=ax, color='navy')
la_highways.plot(ax=ax, color='maroon')
Explanation: Combining the results
Now that we have made a number of queries and joins, let's combine them into a single plot.
To make the plot a bit nicer, we can also load some shapefiles for the coast and land:
End of explanation
<END_TASK>