Geopandas

Geodataframe from a shape

import geopandas as gpd
from tabulate import tabulate
 
myShape = 'C:\\Path\\Of\\My\\Shape.shp'
 
print('\n' + myShape)
 
df = gpd.read_file(myShape, encoding='iso-8859-1')
df['type'] = df['geometry'].astype(str).str.replace(r' .*|\(.*', '', regex=True)
df = df[['id', 'type', 'geometry']]
 
print(tabulate(df.head(10), headers='keys', tablefmt='psql', showindex=True))
print(df.shape[0])
 
MyFieldList = df['id'].drop_duplicates().dropna().sort_values(ascending=False).tolist()
print('\n' + str(MyFieldList))
 
MyGeomTypeList = df['type'].drop_duplicates().dropna().sort_values(ascending=False).tolist()
print('\n' + str(MyGeomTypeList))

Geodataframe from a dataframe using WKT

First, read a dataframe with WKT data in field name wkt for example. Then:

import geopandas as gpd
from shapely import wkt
 
...
 
df['geometry'] = gpd.GeoSeries.from_wkt(df['wkt'])
 
geo_df = gpd.GeoDataFrame(df, geometry='geometry')
geo_df.drop('wkt', axis=1, inplace=True)
 
...

Add CRS/projection

df = df.set_crs('epsg:4326')

Convert CRS/projection

df = df.to_crs(2154)

Display a map

from tabulate import tabulate
import matplotlib.pyplot as plt
 
MyLayer = 'E:\\Path\\to\\shape.shp'
df = gpd.read_file(MyLayer)
 
df.plot()
plt.title('My layer', pad=10, fontsize=10)
plt.show()

Display several geodataframes in one map

fig, ax = plt.subplots()
 
# POLYGONS
df1.plot(
ax=ax,
color='#90EE90',
edgecolor='black',
zorder=2
)
 
# POLYGONS
df2.plot(
ax=ax,
color='#FF7F7F',
edgecolor='black',
zorder=1
)
 
# POINTS
dfPoints.plot(
ax=ax,
facecolor='red',
edgecolor='black',
zorder=3
)
 
# LINES
dfLine.plot(
ax=ax,
edgecolor='purple',
linewidth=10,
zorder=3
)
 
 
plt.title('My map', pad=10, fontsize=10)
plt.show()

Zoom on a layer from a subplot

...
# Extract extent from your GDF
MinX, MinY, MaxX, MaxY = gdf.geometry.total_bounds
print('MinX :', MinX)
print('MinY :', MinY)
print('MaxX :', MaxX)
print('MaxY :', MaxY)
...
# Zoom 
ax.set_xlim(cibleMinX - 1, cibleMaxX + 1)
ax.set_ylim(cibleMinY - 1, cibleMaxY + 1)
...

Display the Coordinate Reference Systems (CRS)

print('Proj : ' + str(df.crs))

Check geometry

df['valid ?'] = df.is_valid
df = df[df['valid ?'] == False]
 
print(tabulate(df.head(5), headers='keys', tablefmt='psql', showindex=False))
print(df.shape[0])

Check projection

if str(df.crs) == 'epsg:2154':
    print(colored('OK: ' + str(df.crs), 'green'))
if str(df.crs) != 'epsg:2154':
    print(colored('Warning: ' + str(df.crs), 'red'))

Extract vertices (all vertex/point) from polygon

Use extract_unique_points() to extract as MultiPoints:

dfMultiPoints = df.copy()
dfMultiPoints['geometry'] = dfMultiPoints['geometry'].extract_unique_points()

Count single geometry

df['Number'] = df['geometry'].apply(lambda geom: len(geom.geoms))

Simplify polygon

import topojson as tp
 
topo = tp.Topology(df, prequantize=False)
 
dfSimple = topo.toposimplify(5).to_gdf()
 
print(tabulate(dfSimple, headers='keys', tablefmt='psql', showindex=False))

Extract segments from polygon

from shapely.geometry import LineString
 
myPolygon = myGdf.geometry.iloc[0]
pointsCoords = list(myPolygon.exterior.coords)
segments = [LineString([pointsCoords[i], pointsCoords[i + 1]]) for i in range(len(pointsCoords) - 1)]
 
dfSegments = gpd.GeoDataFrame(geometry=segments, crs=myGdf.crs)
 
print(tabulate(dfSegments, headers='keys', tablefmt='psql', showindex=True))

Calculate intersection lines on polygons

def calculate_intersection_ratio(line, polygons_union):
    intersection = line.intersection(polygons_union)
    total_length = line.length
    intersected_length = intersection.length if not intersection.is_empty else 0
    return (intersected_length / total_length) * 100 if total_length > 0 else 0
 
polygons_union = dfPolygons.geometry.unary_union
 
dfLines['intersection'] = dfLines.geometry.apply(
    lambda line: calculate_intersection_ratio(line, polygons_union)
)
 
dfLines = dfLines[['Seg Id', 'intersection', 'geometry']]
 
print(tabulate(dfLines, headers='keys', tablefmt='psql', showindex=False))

Calculate line length

dfLines['longueur'] = dfLines.geometry.length

Create an empty shape

import geopandas as gpd
 
myShape = 'C:\\Path\\beautiful shape.shp'
schema = {"geometry": "Polygon", "properties": {"myid": "int"}}
crs = "EPSG:2154"
dfEmpty = gpd.GeoDataFrame(geometry=[])
dfEmpty.to_file(myShape, driver='ESRI Shapefile', schema=schema, crs=crs)

Save a map

import matplotlib.pyplot as plt
from PIL import Image, ImageOps
 
imagePath = myDirectory + '_Carte.png'
plt.savefig(imageChemin, dpi=250)
im = Image.open(imagePath)
bordered = ImageOps.expand(im, border=1, fill=(0, 0, 0))
bordered.save(imagePath)
# im.show()

Hide the axes graduations

ax.xaxis.set_ticks([])
ax.xaxis.set_ticklabels([])
 
ax.yaxis.set_ticks([])
ax.yaxis.set_ticklabels([])

Build a work environment

You are processing geometries in a loop. You want to display sample maps of 6 randomly processed geometries simultaneously.

Start by exporting your maps as an image, and save the paths to your images in a list.

...
imagesTemplate = []
...
imagesTemplate.append(imagePath)
...

Then:

...
from PIL import Image, ImageOps, ImageTk
from tkinter import Tk, Toplevel, Label, PhotoImage
...
# Settings for arranging windows
screen_width = 1920 # Screen width
screen_height = 1020 # Screen height
 
columns = 3 
rows = 2 
 
#Size windows to fit on screen
window_width = screen_width // columns
 
# Adjust
adjustment_factor = 0.905 # Réduire la hauteur des fenêtres de la première ligne
window_height_first_line = int((screen_height // rows) * adjustment_factor) # For 1st line
window_height_rest = screen_height // rows # For 2nd line
 
root = Tk()
root.withdraw() # Cache la fenêtre principale de tkinter
 
# Function to display an image in a tkinter window
def open_image(path, x, y, height):
    top = Toplevel()
    top.geometry(f"{window_width}x{height}+{x}+{y}")
    top.title(path)
 
    img = Image.open(path)
    img = img.resize((window_width, height), Image.LANCZOS)
    img_tk = ImageTk.PhotoImage(img)
 
    label = Label(top, image=img_tk)
    label.image = img_tk
    label.pack()
 
for index, image_path in enumerate(imagesTemplate[0:6]):
    col = index % columns
    row = index // columns
    x = (col * window_width)-8
 
    if row == 0: # 1st line
        y = row * window_height_first_line
        height = window_height_first_line
 
    else: # 2nd line
        y = row * window_height_rest
        y -= 18
        height = window_height_rest
 
    open_image(image_path, x, y, height)
 
root.mainloop()

Or with 2 axes (a map to the left, and some data/text to the right): 

fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [2, 1]})
 
ax[0].set_aspect('equal')
ax[0].set_title('Blablbla' + str(idParcelle) + ' (ID ' + str(id) + ')', pad=10, fontsize=10, weight='bold')
 
dfCible.plot(
ax=ax[0],
color='#0088CC',
edgecolor='black',
zorder=2
)
 
ax[1].axis('off')
ax[1].set_title('Blablabla Xxxx', pad=10, fontsize=10, weight='bold')

Liens ou pièces jointes

France-Departements-Deformation.zip	[France-Departements-Deformation]	335 Ko
simple_countries.zip	[simple_countries]	1880 Ko