2.1 GeoPython - Praktická ukázka¶

Vzorek dat ke stažení zde.

GeoPandas¶

Úvod do práce s datovými typy knihovny GeoPandas najdete v rámci předmětu UZPR.

Vstupní data načteme pomocí metody read_file().

In [15]:

data = gpd.read_file("data/polygons.shp")
data

data = gpd.read_file("data/polygons.shp") data

Out[15]:

	label	note	geometry
0	1	glacier	POLYGON ((630019.737 7747367.105, 630896.053 7...
1	2	debries	POLYGON ((631944.079 7740883.553, 632169.079 7...
2	3	rock	POLYGON ((625330.263 7741517.105, 625780.263 7...

In [16]:

data.set_index("label", inplace=True)
data

data.set_index("label", inplace=True) data

Out[16]:

	note	geometry
label
1	glacier	POLYGON ((630019.737 7747367.105, 630896.053 7...
2	debries	POLYGON ((631944.079 7740883.553, 632169.079 7...
3	rock	POLYGON ((625330.263 7741517.105, 625780.263 7...

In [17]:

data.loc[1]["geometry"]

data.loc[1]["geometry"]

Out[17]:

In [18]:

data.area

data.area

Out[18]:

label
1    3.099765e+06
2    1.513141e+05
3    6.764958e+05
dtype: float64

In [19]:

data.centroid

data.centroid

Out[19]:

label
1      POINT (630375.4 7746334.905)
2    POINT (631925.346 7740645.252)
3    POINT (625216.784 7741106.275)
dtype: geometry

In [20]:

data.boundary

data.boundary

Out[20]:

label
1    LINESTRING (630019.737 7747367.105, 630896.053...
2    LINESTRING (631944.079 7740883.553, 632169.079...
3    LINESTRING (625330.263 7741517.105, 625780.263...
dtype: geometry

In [21]:

data.buffer(100)

data.buffer(100)

Out[21]:

label
1    POLYGON ((629977.831 7747457.901, 629986.163 7...
2    POLYGON ((631897.596 7740972.092, 631906.656 7...
3    POLYGON ((625298.64 7741611.974, 625308.345 77...
dtype: geometry

Modifikace datových objetků¶

Přidáme do datové struktury nové sloupce.

In [22]:

data["area_km2"] = data.area / 1e6
data

data["area_km2"] = data.area / 1e6 data

Out[22]:

	note	geometry	area_km2
label
1	glacier	POLYGON ((630019.737 7747367.105, 630896.053 7...	3.099765
2	debries	POLYGON ((631944.079 7740883.553, 632169.079 7...	0.151314
3	rock	POLYGON ((625330.263 7741517.105, 625780.263 7...	0.676496

In [23]:

data["centroid"] = data.centroid
data

data["centroid"] = data.centroid data

Out[23]:

	note	geometry	area_km2	centroid
label
1	glacier	POLYGON ((630019.737 7747367.105, 630896.053 7...	3.099765	POINT (630375.4 7746334.905)
2	debries	POLYGON ((631944.079 7740883.553, 632169.079 7...	0.151314	POINT (631925.346 7740645.252)
3	rock	POLYGON ((625330.263 7741517.105, 625780.263 7...	0.676496	POINT (625216.784 7741106.275)

In [24]:

first_cen = data["centroid"].iloc[0]
data["distance"] = data["centroid"].distance(first_cen)
data

first_cen = data["centroid"].iloc[0] data["distance"] = data["centroid"].distance(first_cen) data

Out[24]:

	note	geometry	area_km2	centroid	distance
label
1	glacier	POLYGON ((630019.737 7747367.105, 630896.053 7...	3.099765	POINT (630375.4 7746334.905)	0.000000
2	debries	POLYGON ((631944.079 7740883.553, 632169.079 7...	0.151314	POINT (631925.346 7740645.252)	5896.989221
3	rock	POLYGON ((625330.263 7741517.105, 625780.263 7...	0.676496	POINT (625216.784 7741106.275)	7345.058814

In [25]:

data.plot("area_km2", legend=True, figsize=(10,10))

data.plot("area_km2", legend=True, figsize=(10,10))

Out[25]:

<Axes: >

Export¶

Zapsat GeoDataFrame do výstupního souboru lze pomocí funkce to_file():

In [26]:

data[["note", "centroid", "distance"]].to_file(
    "/tmp/polygons_modified.geojson", driver="GeoJSON")

data[["note", "centroid", "distance"]].to_file( "/tmp/polygons_modified.geojson", driver="GeoJSON")

Rasterio¶

Nejprve otevřeme ukázkový datový soubor a vytiskneme jeho metadata:

In [27]:

import rasterio as rio

with rio.open("data/landsat.tif") as ds:
    meta = ds.meta
    print(meta)

import rasterio as rio with rio.open("data/landsat.tif") as ds: meta = ds.meta print(meta)

{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 500, 'height': 500, 'count': 6, 'crs': CRS.from_epsg(32621), 'transform': Affine(30.0, 0.0, 624000.0,
       0.0, -30.0, 7750000.0)}

Rasterizace¶

V dalším kroku provedeme rasterizaci polygonů z minulé části.

In [28]:

from rasterio import features

meta["count"] = 1
polygons_rast_fn = "/tmp/polygons_rast.tif"
# w+: kromě zápisu umožňuje i čtení ze souboru.
with rio.open(polygons_rast_fn, "w+", **meta) as out:
    arr = out.read(1)

    shapes = []
    for geom, value in zip(data.geometry, data.index):
        shapes.append((geom, value))
    # or
    # shapes = ((geom, value) for geom, value in zip(data.geometry, data.index))

    burned = features.rasterize(shapes=shapes, fill=0, out=arr, transform=out.transform)
    out.write_band(1, burned)

from rasterio import features meta["count"] = 1 polygons_rast_fn = "/tmp/polygons_rast.tif" # w+: kromě zápisu umožňuje i čtení ze souboru. with rio.open(polygons_rast_fn, "w+", **meta) as out: arr = out.read(1) shapes = [] for geom, value in zip(data.geometry, data.index): shapes.append((geom, value)) # or # shapes = ((geom, value) for geom, value in zip(data.geometry, data.index)) burned = features.rasterize(shapes=shapes, fill=0, out=arr, transform=out.transform) out.write_band(1, burned)

Výsledek rasterizace zkontrolujeme:

In [29]:

from pathlib import Path

Path(polygons_rast_fn).exists()

from pathlib import Path Path(polygons_rast_fn).exists()

Out[29]:

True

In [30]:

with rio.open(polygons_rast_fn) as ds:
    poly_arr = ds.read(1)
    print(poly_arr.shape)

with rio.open(polygons_rast_fn) as ds: poly_arr = ds.read(1) print(poly_arr.shape)

(500, 500)

In [31]:

import numpy as np
np.unique(poly_arr, return_counts=True)

import numpy as np np.unique(poly_arr, return_counts=True)

Out[31]:

(array([0, 1, 2, 3], dtype=int16), array([245630,   3446,    167,    757]))

In [32]:

import matplotlib.pyplot as plt

plt.imshow(poly_arr)
plt.colorbar()

import matplotlib.pyplot as plt plt.imshow(poly_arr) plt.colorbar()

Out[32]:

<matplotlib.colorbar.Colorbar at 0x7eeebc993640>

Operace překrytí¶

Nejprve načteme vstupní vrstvu výřezu snímku Landsat:

In [33]:

import rasterio as rio

with rio.open("data/landsat.tif") as ds:
    print(ds.count)
    print(ds.indexes)
    image_arr = np.dstack([ds.read(b) for b in ds.indexes])

import rasterio as rio with rio.open("data/landsat.tif") as ds: print(ds.count) print(ds.indexes) image_arr = np.dstack([ds.read(b) for b in ds.indexes])

6
(1, 2, 3, 4, 5, 6)

In [34]:

image_arr.shape, poly_arr.shape

image_arr.shape, poly_arr.shape

Out[34]:

((500, 500, 6), (500, 500))

Data vykreslíme:

In [35]:

plt.subplot(1, 3, 1)
plt.imshow(image_arr[:, :, 2], cmap=plt.cm.Greys_r)
plt.title('Red band')

plt.subplot(1, 3, 2)
plt.imshow(image_arr[:, :, 5], cmap=plt.cm.Greys_r)
plt.title('SWIR band')

plt.subplot(1, 3, 3)
plt.imshow(poly_arr, cmap=plt.cm.Spectral, interpolation='nearest')
plt.title('Training Data')

plt.show()

plt.subplot(1, 3, 1) plt.imshow(image_arr[:, :, 2], cmap=plt.cm.Greys_r) plt.title('Red band') plt.subplot(1, 3, 2) plt.imshow(image_arr[:, :, 5], cmap=plt.cm.Greys_r) plt.title('SWIR band') plt.subplot(1, 3, 3) plt.imshow(poly_arr, cmap=plt.cm.Spectral, interpolation='nearest') plt.title('Training Data') plt.show()

Na závěr si vyzkoušíme převedení hodnot do jednorozměrného pole. Tuto operaci využijeme v následujících lekcích.

In [36]:

poly_arr_flat = poly_arr[poly_arr > 0]
poly_arr_flat.shape

poly_arr_flat = poly_arr[poly_arr > 0] poly_arr_flat.shape

Out[36]:

(4370,)

In [37]:

np.unique(poly_arr_flat, return_counts=True)

np.unique(poly_arr_flat, return_counts=True)

Out[37]:

(array([1, 2, 3], dtype=int16), array([3446,  167,  757]))

Porovnejte s:

In [38]:

np.unique(poly_arr, return_counts=True)

np.unique(poly_arr, return_counts=True)

Out[38]:

(array([0, 1, 2, 3], dtype=int16), array([245630,   3446,    167,    757]))

Takto vytvořené pole lze použít jako masku pro obrazová data.

In [39]:

image_arr[poly_arr > 0, :].shape

image_arr[poly_arr > 0, :].shape

Out[39]:

(4370, 6)

Na závěr připravíme data pro další kroky strojového účení.

In [41]:

x = image_arr[poly_arr > 0]
y = poly_arr[poly_arr > 0]
x.shape, y.shape

x = image_arr[poly_arr > 0] y = poly_arr[poly_arr > 0] x.shape, y.shape

Out[41]:

((4370, 6), (4370,))