Das LGLN stellt zahlreiche offene Geodaten bereit. Diese Jupyter-Notebook-Serie demonstriert den Zugriff mit Python.
https://github.com/BostelmannLGLN/LGLN-OpenData-Notebooks/blob/main/requirements.txt
#!pip install -r requirements.txt
from pathlib import Path
import rasterio
from rasterio.windows import Window
from rasterio.transform import Affine
from rasterio.plot import reshape_as_raster, reshape_as_image, show
from skimage import img_as_ubyte
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
from random import randrange
from IPython.display import Video
import animation
DOPs = gpd.read_file('https://single-datasets.s3.eu-de.cloud-object-storage.appdomain.cloud/pro-download-indices/dop/lgln-opengeodata-dop20.geojson')
DOPs
DOPs.plot(edgecolor='black', figsize=(15, 15))
rgb_image_url = DOPs.iloc[randrange(len(DOPs))].rgb
rgb_image_url
links = 342000
rechts = 384000
unten = 5826000
oben = 5838000
DOPs.cx[links:rechts, unten:oben].plot(edgecolor='black')
multi_DOPs = DOPs[DOPs['tile_id'].duplicated()]
#multi_DOPs.to_file('multi_DOPs.geojson', driver='GeoJSON') # In Datei schreiben
base = DOPs.plot(figsize=(15, 15))
multi_DOPs.plot(ax=base, color='orange')
multi_DOPs
double_ids = multi_DOPs.tile_id
tile_id = double_ids.iloc[randrange(len(multi_DOPs))]
tile_id
tile_id = '324865862' # Syke
rgb_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgb.to_list()
rgb_url_list
%%time
with rasterio.open(rgb_url_list[0]) as cog: # Cloud Optimized GeoTiff
print(cog.profile)
rgb_data = cog.read()
plt.figure(figsize=(15, 15))
ax = show(rgb_data)
%%time
with rasterio.open(rgb_url_list[1]) as cog: # Cloud Optimized GeoTiff
print(cog.profile)
rgb_data = cog.read(out_shape=(3, 1000, 1000))
#(Kanäle, Höhe, Breite))
plt.figure(figsize=(15, 15))
ax = show(rgb_data)
pixel_von_oben = 6000
pixel_von_links = 4000
hoch = 1000
breit = 1000
images = []
dates = []
for url in rgb_url_list:
with rasterio.open(url) as cog:
profile = cog.profile
profile.update({"driver": "GTiff",
"height": hoch,
"width": breit})
rgb_data = cog.read(out_shape=(3, hoch, breit),
window = Window(pixel_von_links,
pixel_von_oben,
breit,
hoch))
#rasterio.windows.Window(col_off, row_off, width, height)
# Anzeigen
plt.figure(figsize=(6, 6))
ax = show(rgb_data)
# Umbennen
date = url.split('/')[4].replace("-", "")
dates.append(date)
file_name = f'dop_{date}_rgb.tif'
images.append(file_name)
# Schreiben
with rasterio.open(file_name, "w", **profile) as output:
output.write(rgb_data)
print(dates, images)
out_path = animation.create_animation(images, dates,
name=tile_id,
type='mp4',
duration=1,
add_ban=True,
add_bar=False,
add_name=True,
center_text=f'{int(breit*cog.res[0])} m x {int(hoch*cog.res[1])} m',
logo_path=(Path('./logos/logo_left.png'),
Path('./logos/no_logo_right.png')))
Video(out_path, html_attributes='loop autoplay')
rgbi_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgbi.to_list()
rgbi_url = rgbi_url_list[1]
date = rgbi_url.split('/')[4].replace("-", "")
date
cog = rasterio.open(rgbi_url)
rgbi = cog.read(out_shape=(4, hoch, breit),
window = Window(pixel_von_links,
pixel_von_oben,
breit,
hoch)).astype('float32')
red = rgbi[0]
green = rgbi[1]
blue = rgbi[2]
nir = rgbi[3]
plt.figure(figsize=(15, 15))
ax = show(nir, cmap='gray')
np.seterr(divide='ignore', invalid='ignore')
ndvi = np.where(
(nir+red)==0,
0,
(nir-red)/(nir+red))
ndvi
ndvi.min(), ndvi.max()
plt.figure(figsize=(15, 15))
ax = show(ndvi, cmap='RdYlGn')
out_profile = cog.profile
t = out_profile['transform']
t = Affine(t[0],
t[1],
t[2]+pixel_von_links*t[0],
t[3],
t[4],
t[5]+pixel_von_oben*t[4])
out_profile.update(
{"dtype": 'float32',
"count": 1,
"height": hoch,
"width": breit,
"transform": t})
out_profile
with rasterio.open(f'dop_{date}_ndvi.tif', "w", **out_profile) as output:
output.write(ndvi, 1)
def get_colors(inp, colormap, vmin=None, vmax=None):
norm = plt.Normalize(vmin, vmax)
return colormap(norm(inp))[0]
def plot(in_path, out_path, driver='GTiff', cm=plt.cm.viridis , min_value=0, max_value=255):
in_image = rasterio.open(in_path)
in_profile = in_image.profile
in_data = in_image.read()
out_profile = in_profile
out_profile.update(
{"driver": driver,
"count": 3,
"dtype": 'uint8'})
out_data = get_colors(in_data, cm, min_value, max_value)
out_data = img_as_ubyte(out_data[...,:-1])
out_data = reshape_as_raster(out_data)
with rasterio.open(out_path, "w", **out_profile) as output:
output.write(out_data)
plot(f'dop_{date}_ndvi.tif', f'dop_{date}_ndvi_plot.tif', 'GTIFF', plt.cm.RdYlGn, ndvi.min(), ndvi.max())
plot(f'dop_{date}_ndvi.tif', f'dop_{date}_ndvi_plot.png', 'PNG', plt.cm.RdYlGn, ndvi.min(), ndvi.max())
rgbi_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgbi.to_list()
def ndvi(rgbi_url, pixel_von_oben=0, pixel_von_links=0, breit=1000, hoch=1000):
date = rgbi_url.split('/')[4].replace("-", "")
cog = rasterio.open(rgbi_url)
rgbi = cog.read(out_shape=(4, hoch, breit),
window = Window(pixel_von_links,
pixel_von_oben,
breit,
hoch)).astype('float32')
red = rgbi[0]
nir = rgbi[3]
np.seterr(divide='ignore', invalid='ignore')
ndvi = np.where(
(nir+red)==0,
0,
(nir-red)/(nir+red))
return ndvi
ndvi_1 = ndvi(rgbi_url_list[0], pixel_von_oben=pixel_von_oben, pixel_von_links=pixel_von_links, breit=breit, hoch=hoch)
# Bei Fehlermeldung nochmal probieren
ndvi_2 = ndvi(rgbi_url_list[1], pixel_von_oben=pixel_von_oben, pixel_von_links=pixel_von_links, breit=breit, hoch=hoch)
# Bei Fehlermeldung nochmal probieren
ndvi_diff = ndvi_2 - ndvi_1
plt.figure(figsize=(16, 16))
ax = show(ndvi_diff, cmap='BrBG')
with rasterio.open(f'dop_{date}_ndvi_diff.tif', "w", **out_profile) as output:
output.write(ndvi_diff, 1)
plot(f'dop_{date}_ndvi_diff.tif', f'dop_{date}_ndvi_diff_plot.tif',
'GTiff', plt.cm.BrBG, ndvi_diff.min(), ndvi_diff.max())