In [136]:
import fiona
import geopandas as gpd
from osgeo import gdal
from osgeo import ogr
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
Fiona is a minimalist python package for reading (and writing) vector data in python. Fiona provides python objects (e.g. a dictionary for each record) to geospatial data in various formats.
Reading in vector files with Fiona is - good practice - done using the with environment:
In [55]:
with fiona.open('../data/deelbekkens/Deelbekken.shp') as deelbekkens:
feature = next(iter(deelbekkens)) # Just one checking the first
print("Bekken: ", feature['properties']['BEKNAAM'])
print("Vectortype: ", feature['geometry']['type'])
print(feature['geometry']['coordinates'][0][0])
This is actually equivalent to:
In [3]:
deelbekkens = fiona.open('../data/deelbekkens/Deelbekken.shp')
feature = next(iter(deelbekkens))
print("Bekken: ", feature['properties']['BEKNAAM'])
print("Stroomgebied: ", feature['properties']['STRMGEB'])
print("Vectortype: ", feature['geometry']['type'])
deelbekkens.close()
But people tend to forget to use the close(!)
The projection information is interpreted by Fiona as well:
In [33]:
with fiona.open('../data/deelbekkens/Deelbekken.shp') as deelbekkens:
#print(deelbekkens.crs_wkt)
print(deelbekkens.crs)
print(deelbekkens.crs_wkt)
Actually, **GeoPandas relies on Fiona** and translates this information into a tabular information, compatible to Pandas:
In [5]:
df_deelbekkens = gpd.read_file('../data/deelbekkens/Deelbekken.shp')
df_deelbekkens.head()
Out[5]:
The capabilities to read file with Fiona is dependent on the available drivers on your computer. You can get an overview by asking fiona for the supported drivers:
In [6]:
fiona.supported_drivers
Out[6]:
It provides also information about what you can do with the different types of data:
r: readinga: appending (adding elements)w: writing
In [7]:
with fiona.open('../data/EUgrid10.geojson') as eugrid:
feature = next(iter(eugrid)) # Just one checking the first
print("cellcode: ", feature['properties']['CellCode'])
print("Vectortype: ", feature['geometry']['type'])
In [8]:
gpd.read_file('../data/EUgrid10.geojson').head(1)
Out[8]:
So, what about GDAL?
In [9]:
driver = ogr.GetDriverByName('ESRI Shapefile')
deelbekken_gdal = driver.Open('../data/deelbekkens/Deelbekken.shp')
deelbekken_gdal.GetLayerCount()
Out[9]:
In [10]:
deelbekken_layer = deelbekken_gdal.GetLayer()
print("The number of features of the deelbekken shapefile is: ", deelbekken_layer.GetFeatureCount())
That's why they call it not pythonic...
However, GDAL is powerfull for other activities, we'll come back to that later. Moreover, it supports mostly more file formats as fiona (GeoPandas) provide:
In [11]:
!ogr2ogr --formats
For example, the interaction with an Esri File GeoDataBase, is provided in this example...
Shapely is a Python library for geometric operations using the **GEOS** library.
Shapely can perform:
In [12]:
from shapely.geometry import LineString, Polygon
In [13]:
line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
line
Out[13]:
In [14]:
dilated = line.buffer(0.5, cap_style=3)
dilated
Out[14]:
Getting back the coorindates as x, y arrays:
In [15]:
polygon = Polygon(LineString([(0.2, 0.2), (0.2, 1), (1, 1), (1, 0.2)]))
polygon
Out[15]:
In [16]:
line.coords.xy
Out[16]:
Binary predicates to check if the line is contained by the polygon
In [17]:
line.within(polygon)
Out[17]:
In [29]:
polygon.intersects(line)
Out[29]:
Actually, **GeoPandas relies on shapely** and uses these geometric properties and operations
In [19]:
deelbekkens = gpd.read_file("../data/deelbekkens/Deelbekken.shp")
In [22]:
demer = deelbekkens[deelbekkens["BEKNAAM"] == "Demerbekken"]
In [28]:
demer.plot()
demer.convex_hull.plot()
Out[28]:
Pyproj handles the Spatial Reference Systems (SRS) transformations. Actually, the only items to remember are:
pyproj.Proj: define a SRSpyproj.transform: execute a transformation between two defined SRS
In [3]:
from pyproj import Proj, transform
The easiest way is mostly just using the EPSG code:
Lambert 72:
In [4]:
srs_lambert = Proj(init='epsg:31370')
In [5]:
srs_lambert.srs
Out[5]:
WGS84:
In [6]:
srs_wgs84 = Proj(init='epsg:4326')
In [72]:
srs_wgs84.is_latlong()
Out[72]:
Other methods to define the SRS are also available:
In [65]:
utm32 = Proj(proj="utm", zone="32")
The transform function takes the two defined SRS objects, together with an X/Y coordinate combination
In [79]:
transform(srs_lambert, srs_wgs84, 98710.32800000161, 162573.7030000016)
Out[79]:
Actually, **GeoPandas relies on pyproj** to perform SRS transformations using the EPSG code (but calls it `coordinate reference system (CRS)`...)
In [83]:
deelbekkens = gpd.read_file("../data/deelbekkens/Deelbekken.shp")
deelbekkens.crs
Out[83]:
In [170]:
deelbekkens.to_crs(epsg="4326").head()
Out[170]:
Actually, the function is useful in general to transform the X/Y combinations in any coordinate table. Consider the following example of species data with the coordinates in WGS84.
In [106]:
inv_data = pd.read_excel("../data/invasive_extract.xlsx")
In [107]:
inv_data.head()
Out[107]:
Conversion to Lambert 72 is supported by the pyproj Package as follows. As I want to support future usage of this functionality as well, I'll write my own custom function:
In [125]:
def transform_dfrow_projection(row, x_name, y_name, srs_from, srs_to):
"""
Converts the x and y coordinates of a given DataFrame row into a Series of the
longitude and latitude.
"""
from pyproj import transform
return pd.Series(transform(srs_from, srs_to, row[x_name], row[y_name]))
For a single row, this looks like:
In [126]:
transform_dfrow_projection(inv_data.loc[0], "decimalLongitude", "decimalLatitude",
srs_wgs84, srs_lambert)
Out[126]:
Applying this to the entire table:
In [127]:
inv_data[["x", "y"]] = inv_data.apply(transform_dfrow_projection, axis=1,
args=("decimalLongitude", "decimalLatitude",
srs_wgs84, srs_lambert))
This feels too complicated?!?Understandable ,but let's just once take the effort to compile this into a reusable function...
In [130]:
def add_trf_coordinates(filename, x_name="decimalLongitude", y_name="decimalLatitude",
srs_from='epsg:4326', srs_to='epsg:31370',
new_x="X", new_y="Y"):
"""Add new X/Y columns with transformed SRS to a given Excel dataset with X/Y columns
Parameters
----------
filename : str
Path to excel file to read dat from
x_name : str
Column containing the x-coordinates (default decimalLongitude)
y_name : str
Column containing the y-coordinates (default decimalLatitude)
srs_from: str
epsg:n code for the current SRS of the coordinates
srs_to: str
epsg:n code for the SRS of the added coordinates
new_x: str
name of the added x-coordinate column
new_y: str
name of the added y-coordinate column
returns
-------
df with 2 additional columns
"""
df = pd.read_excel("../data/invasive_extract.xlsx")
srs_from = Proj(init=srs_from)
srs_to = Proj(init=srs_to)
df[[new_x, new_y]] = df.apply(transform_dfrow_projection, axis=1,
args=(x_name, y_name, srs_from, srs_to))
return df
In [134]:
df_enriched = add_trf_coordinates("../data/invasive_extract.xlsx")
df_enriched.head()
Out[134]:
Remark that GDAL also has a utility gdalsrsinfo to get/check the SRS information of GIS file:
In [132]:
! gdalsrsinfo -o proj4 ../data/EUgrid10.geojson #try also wiht -o all
Or to execute the transformation of a dataset together with the transformation of the data type:
In [1]:
!ogr2ogr -t_srs "EPSG:4326" -f "GeoJSON" "../scratch/deelbekken.geojson" "../data/deelbekkens/Deelbekken.shp"
(More on this in the next notebook)
Actually, we got a whole sequence of data types (objects) to end up with GeoDataFrames:
Our df_enriched is a Pandas DataFrame, so we can convert to a GeoDataFrame by converting the coordinates to a set of POINT geometries:
In [217]:
from shapely.geometry import Point
geometry = [Point(xy) for xy in zip(df_enriched["X"], df_enriched["Y"])]
crs = {'init': 'epsg:31370'}
geo_df = gpd.GeoDataFrame(df_enriched, crs=deelbekkens.crs, geometry=geometry)
In [197]:
import mplleaflet
In [214]:
fig, ax = plt.subplots()
geo_df.to_crs(epsg="4326").plot(markersize=10, ax=ax) # on the fly conversion to WGS84
mplleaflet.display()
Out[214]:
In [208]:
gpd.sjoin(geo_df, deelbekkens[["BEKNAAM", "UIDN", "geometry"]], how='left', op="within")
Out[208]:
In [215]:
fig, ax = plt.subplots()
geo_df.to_crs(epsg="4326").plot(markersize=10, ax=ax) # on the fly conversion to WGS84
geo_df.buffer(500).to_crs(epsg="4326").plot(ax=ax) # on the fly conversion to WGS84
mplleaflet.display()
Out[215]:
For those interested, also check the following example from Joris Vandenbossche, explaining how to derive the distance to the nearest green area for all adresses in the centre of Antwerp: