2  Vector layers

2.1 Loading packages

First, we import the two packages we will be working with:

• pandas—for working with tables, and
• geopandas—For working with vector layers,

as follows:

import pandas as pd
import geopandas as gpd

2.2 Table to point layer

In the first example, we are going to read an Excel file with x/y coordinates, and convert it to a (point) vector layer. The file bycode2019-2.xlsx contains town locations in Israel in 2019, downloaded from data.gov.il.

To read an Excel (.xlsx) file, we use the pd.read_excel function. The result is a data structure called DataFrame, used to represent a table in the pandas package:

dat = pd.read_excel('data/bycode2019-2.xlsx')
dat

	שם יישוב	סמל יישוב	תעתיק	...	שנה	שם יישוב באנגלית	אשכול רשויות מקומיות
0	אבו ג'ווייעד (שבט)	967	ABU JUWEI'ID	...	2019	Abu Juway'ad	NaN
1	אבו גוש	472	ABU GHOSH	...	2019	Abu Ghosh	NaN
2	אבו סנאן	473	ABU SINAN	...	2019	Abu Sinan	NaN
...	...	...	...	...	...	...	...
1480	תראבין א-צאנע (שבט)	970	TARABIN AS-SANI	...	2019	Tarabin as-Sani'	NaN
1481	תרבין א-צאנע (יישוב)*	1346	TARABIN AS-SANI	...	2019	Tarabin As-Sani	610.0
1482	תרום	778	TARUM	...	2019	Tarum	NaN

1483 rows × 23 columns

The column named “coordinates” contains the town coordinates:

dat['קואורדינטות']

0       2.040057e+09
1       2.105263e+09
2       2.160776e+09
            ...     
1480    1.830056e+09
1481    1.752658e+09
1482    1.983663e+09
Name: קואורדינטות, Length: 1483, dtype: float64

The following code section:

• filters out the rows with missing coordinate values, since these cannot be tranalated to point geometries, and
• splits the coordinate string to X and Y values, according to the instructions.

dat = dat[dat['קואורדינטות'].notna()].copy().reset_index(drop=True)
dat['x'] = dat['קואורדינטות'].astype(str).str.slice(0,5).astype(float) * 10
dat['y'] = dat['קואורדינטות'].astype(str).str.slice(5,10).astype(float) * 10

As a result, we now have x and y numeric columns with the ITM coordinates of each town:

dat

	שם יישוב	סמל יישוב	תעתיק	...	אשכול רשויות מקומיות	x	y
0	אבו ג'ווייעד (שבט)	967	ABU JUWEI'ID	...	NaN	204000.0	571000.0
1	אבו גוש	472	ABU GHOSH	...	NaN	210520.0	634810.0
2	אבו סנאן	473	ABU SINAN	...	NaN	216070.0	762840.0
...	...	...	...	...	...	...	...
1448	תראבין א-צאנע (שבט)	970	TARABIN AS-SANI	...	NaN	183000.0	564000.0
1449	תרבין א-צאנע (יישוב)*	1346	TARABIN AS-SANI	...	610.0	175260.0	583690.0
1450	תרום	778	TARUM	...	NaN	198360.0	632270.0

1451 rows × 25 columns

The two columns, along with the CRS definition (EPSG:2039 for ITM) can be transformed to a GeoSeries object. A GeoSeries is a sequence of geometries, along with the CRS definition, also known as the “geometry column” when it is part of a vector layer:

geom = gpd.points_from_xy(dat['x'], dat['y'], crs=2039)
geom = gpd.GeoSeries(geom)
geom

0       POINT (204000.000 571000.000)
1       POINT (210520.000 634810.000)
2       POINT (216070.000 762840.000)
                    ...              
1448    POINT (183000.000 564000.000)
1449    POINT (175260.000 583690.000)
1450    POINT (198360.000 632270.000)
Length: 1451, dtype: geometry

One of the first things we might want to do with a GeoSeries is to .plot it, to examine what it looks like:

geom.plot();

Combining a GeoSeries with a corresponding table yields GeoDataFrame object, a data structure representing a vector layer, where:

• The GeoSeries is the geometric/spatial part
• The DataFrame is the non-spatial/attributes part

For example, here we combine the GeoSeries named geom, with a DataFrame named dat, resulting in a GeoDataFrame named pnt:

pnt = gpd.GeoDataFrame(data=dat, geometry=geom)
pnt

	שם יישוב	סמל יישוב	תעתיק	...	x	y	geometry
0	אבו ג'ווייעד (שבט)	967	ABU JUWEI'ID	...	204000.0	571000.0	POINT (204000.000 571000.000)
1	אבו גוש	472	ABU GHOSH	...	210520.0	634810.0	POINT (210520.000 634810.000)
2	אבו סנאן	473	ABU SINAN	...	216070.0	762840.0	POINT (216070.000 762840.000)
...	...	...	...	...	...	...	...
1448	תראבין א-צאנע (שבט)	970	TARABIN AS-SANI	...	183000.0	564000.0	POINT (183000.000 564000.000)
1449	תרבין א-צאנע (יישוב)*	1346	TARABIN AS-SANI	...	175260.0	583690.0	POINT (175260.000 583690.000)
1450	תרום	778	TARUM	...	198360.0	632270.0	POINT (198360.000 632270.000)

1451 rows × 26 columns

Figure 2.1 summarizes the workflow of constructing a vector layer:

• Individual geometries—shapely objects (package shapely)
• Geometry column—GeoSeries objects
• Vector layer—GeoDataFrame object

Figure 2.1: Creating a GeoDataFrame from scratch

Figure 2.2 shows the hierarchical structure of GeoDataFrame, using another hypothetical vector layer with two points (London and Paris):

• The entire vector layer is a GeoDataFrame object
• The geometric part, i.e., the geometry column, is a GeoSeries object
• Each cell in the geometry column, i.e., an individual geometry, is a shapely object (package shapely)

Figure 2.2: Structure of a GeoDataFrame

2.3 Reading from file

Another common workflow with vector layers is to import an existing one, from a file. In the next example, we import a Shapefile named TAZ_NORTH_POPDENS_2016.shp. This is a polygonal layer with population density (POP_DENSE) extimates for Northern Israel. The Shapefile was downloaded from data.gov.il.

To import a vector layer from a file, we use the gpd.read_file function:

pol = gpd.read_file('data/TAZ_NORTH_POPDENS_2016.shp')
pol

	OBJECTID	POP_DENSE	SHAPE_Leng	...	POP2016	TAZ_AREA	geometry
0	1	1425.666870	13215.852473	...	9852	6.910450e+06	POLYGON ((217434.341 73558...
1	2	154.888062	21350.948435	...	2748	1.774185e+07	POLYGON ((215109.748 73483...
2	3	115.392609	20113.127398	...	2064	1.788676e+07	POLYGON ((212411.866 74531...
...	...	...	...	...	...	...	...
778	779	10459.399414	4026.527360	...	6828	6.528099e+05	POLYGON ((227634.472 72406...
779	780	1758.489014	3316.818207	...	1017	5.783374e+05	POLYGON ((226951.788 72387...
780	781	518.559326	20705.740930	...	5568	1.073744e+07	POLYGON ((226132.518 72486...

781 rows × 7 columns

We can see that the geometry type, in this case, is POLYGON.

2.4 Plotting

By default, plotting a GeoDataFrame using .plot shows an image of the geometric part:

pol.plot();

We can add symbology and legend, according to an attribute, using the column and legend=True arguments. For example, here we display the population density estimates (POP_DENSE attributes):

pol.plot(column='POP_DENSE', legend=True);

To modify polygon outline style and the color palette, we can use the following additional arguments:

• edgecolor='black'—Black lines
• linewidth=0.1—Reduced line width
• cmap='Reds—The Reds color palette (from colorbrewer)

pol.plot(
    column='POP_DENSE', 
    legend=True,
    edgecolor='black', 
    linewidth=0.1, 
    cmap='Reds'
);

To display more than one layer in the same plot, we need to:

• store the first plot in a variable (e.g., base), and
• pass it to as the ax argument of any subsequent plot(s) (e.g., ax=base).

For example:

base = pol.plot(color='red', edgecolor='none')
pnt.plot(ax=base, markersize=0.5, color='black');

Another useful method to interactively examine the layer(s) is .explore. This creates an interactive map, with a background basemap for context. The parameters of .explore are mostly similar to .plot, with some differences (such as style_kwds for style settings):

pol.explore(
    column='POP_DENSE', 
    legend=True, 
    cmap='Reds',
    style_kwds={'color': 'black', 'weight': 1, 'fillOpacity': 0.5}
)

Make this Notebook Trusted to load map: File -> Trust Notebook

For more information about mapping with geopandas, see the Mapping and Plotting Tools tutorial.

2.5 Geoprocessing

geopandas provides the standard geoprocessing operators, using shapely and pyproj under the hood (which in turn are interfaces to the GEOS and PROJ software, respectively). For example:

• CRS and reprojection—Transforming a given layer from one CRS to another
• Numeric calculations—Calculating numeric geometry properties, such as length, area, and distance
• New geometries—Creating new geometries, such as calculating buffers, or area of intersection
• Geometric relations—Evaluating the relation between layers, such as whether two geometries intersect
• Spatial join—Joining attributes from one layer to another, based on spatial relations

For example, the .area property gives the area sizes of the geometries in CRS units (\(m^2\)):

pol.area

0      6.910450e+06
1      1.774185e+07
2      1.788676e+07
           ...     
778    6.528099e+05
779    5.783374e+05
780    1.073744e+07
Length: 781, dtype: float64

As a more complicated example, let us calculate the distances between each polygon in pol and one point in pnt (Haifa). First, we subset the pnt feature (Haifa) we are interested in, through “selection by attributes”. The GeoDataFrame named haifa has one feature:

haifa = pnt[pnt['שם יישוב'] == 'חיפה']
haifa

	שם יישוב	סמל יישוב	תעתיק	...	x	y	geometry
544	חיפה	4000	HAIFA	...	201120.0	745440.0	POINT (201120.000 745440.000)

1 rows × 26 columns

From that feature, we can extract the first (and only) “geometry”. This is a shapely object, which by default is plotted:

haifa['geometry'].iloc[0]

Now that we have:

• a vector layer (pol), and
• an individual geometry (haifa['geometry'].iloc[0]),

we can use the .distance method to calculate the distances from all polygons to the point. The result is a numeric Series, which we can immediately “insert” as a new attribute named "dist" in pol:

pol['dist'] = pol.distance(haifa['geometry'].iloc[0])
pol

	OBJECTID	POP_DENSE	SHAPE_Leng	...	TAZ_AREA	geometry	dist
0	1	1425.666870	13215.852473	...	6.910450e+06	POLYGON ((217434.341 73558...	17556.437003
1	2	154.888062	21350.948435	...	1.774185e+07	POLYGON ((215109.748 73483...	14726.885602
2	3	115.392609	20113.127398	...	1.788676e+07	POLYGON ((212411.866 74531...	6678.954051
...	...	...	...	...	...	...	...
778	779	10459.399414	4026.527360	...	6.528099e+05	POLYGON ((227634.472 72406...	33459.684072
779	780	1758.489014	3316.818207	...	5.783374e+05	POLYGON ((226951.788 72387...	33493.912906
780	781	518.559326	20705.740930	...	1.073744e+07	POLYGON ((226132.518 72486...	31011.148069

781 rows × 8 columns

Plotting the resulting dist attribute shows the distribution of distances, which are between \(0\) and \(80,000\) \(m\) (i.e., \(80\) \(km\)):

base = pol.plot(column='dist', legend=True, cmap='Spectral')
haifa.plot(ax=base, color='blue', edgecolor='black', markersize=80);

2.6 Writing to file

To export a GeoDataFrame to a vector layer file, we use the .to_file method. The file format is automatically determined according to the chosen file extension, such as:

• .shp—Shapefile
• .gpkg—GeoPackage
• .geojson—GeoJSON

For example, here is how we can export the pol layer to a GeoPackage file named pol.gpkg:

pol.to_file('pol.gpkg')

Figure 2.3 shows the expored file when opened in QGIS.

Figure 2.3: Exported GeoPackage (.gpkg) file viewed in QGIS

2.7 More information

See the Introduction to GeoPandas tutorial (Figure 2.4), and other sections in the geopandas documentation, for more information about geopandas.

Figure 2.4: The Introduction to GeoPandas tutorial