Komplexní úloha - druhá část¶
Postup:
- Stáhnout DMR5G Geoportál
- Vytvořit DMR v rastrové reprezentaci pomocí PDAL (https://pdal.io/en/2.6.0/pipeline.html#dtm)
- Spočítat rozdíl vodorovné a skloněné plochy
In [ ]:
Copied!
import pdal
import pdal
In [ ]:
Copied!
json = """
[
"/tmp/LITC67.laz",
{
"type": "filters.smrf",
"window": 33,
"slope": 1.0,
"threshold": 0.15,
"cell": 1.0
},
{
"type": "filters.range",
"limits": "Classification[2:2]"
},
{
"type":"writers.gdal",
"filename": "/tmp/LITC67.tif",
"output_type": "min",
"gdaldriver": "GTiff",
"window_size": 3,
"resolution": 1.0
}
]
"""
json = """
[
"/tmp/LITC67.laz",
{
"type": "filters.smrf",
"window": 33,
"slope": 1.0,
"threshold": 0.15,
"cell": 1.0
},
{
"type": "filters.range",
"limits": "Classification[2:2]"
},
{
"type":"writers.gdal",
"filename": "/tmp/LITC67.tif",
"output_type": "min",
"gdaldriver": "GTiff",
"window_size": 3,
"resolution": 1.0
}
]
"""
In [ ]:
Copied!
p = pdal.Pipeline(json)
p.execute()
p = pdal.Pipeline(json)
p.execute()
In [ ]:
Copied!
import fiona
from shapely.geometry import shape
import fiona
from shapely.geometry import shape
In [ ]:
Copied!
# potrebne importy
import rasterio
import rasterio.mask
import numpy as np
import math
# potrebne importy
import rasterio
import rasterio.mask
import numpy as np
import math
In [ ]:
Copied!
# promenne pouzite ve skriptu
parcely_poly = "/tmp/parcely.gml"
# promenne pouzite ve skriptu
parcely_poly = "/tmp/parcely.gml"
In [ ]:
Copied!
schema = {
"geometry": "Polygon",
"properties": {
"cislo": "str:25",
"vymera": "float",
"rozdil": "float"
}
}
schema = {
"geometry": "Polygon",
"properties": {
"cislo": "str:25",
"vymera": "float",
"rozdil": "float"
}
}
In [ ]:
Copied!
tolerance = 0.000001
tolerance = 0.000001
In [ ]:
Copied!
# funkce pro orez rasteru mnohouhelnikem
def clip_raster(raster, vector):
image, _ = rasterio.mask.mask(raster, [vector], crop=True)
data = image.astype(float)[0]
data[data == -9999] = np.nan
return data
# funkce pro orez rasteru mnohouhelnikem
def clip_raster(raster, vector):
image, _ = rasterio.mask.mask(raster, [vector], crop=True)
data = image.astype(float)[0]
data[data == -9999] = np.nan
return data
In [ ]:
Copied!
# funkce pro vypocet prumerneho sklonu
def slope_avg(dem):
dy, dx = np.gradient(dem)
slope = np.sqrt(dx ** 2 + dy ** 2)
return np.nanmean(slope)
# funkce pro vypocet prumerneho sklonu
def slope_avg(dem):
dy, dx = np.gradient(dem)
slope = np.sqrt(dx ** 2 + dy ** 2)
return np.nanmean(slope)
In [ ]:
Copied!
# kolik parcel ma sikmou plochu rozdilnou od vodorovne o vice nez 0.01
with fiona.open("parcely.shp", 'r', driver="ESRI Shapefile",
crs="epsg:5514", schema=schema) as ds:
with rasterio.open('/tmp/LITC67.tif', 'r') as litc:
count = 0
count_dif = 0
for f in ds:
count += 1
try:
a = clip_raster(litc, shape(f['geometry']))
prumerne_prevyseni = slope_avg(a)
sikma_plocha = math.sqrt(
prumerne_prevyseni ** 2 + shape(f["geometry"]).area ** 2
)
rozdil = abs(sikma_plocha - shape(f["geometry"]).area)
if rozdil > tolerance:
count_dif += 1
except ValueError:
pass
print("-" * 80)
print("Počet parcel: {}".format(count))
print("Počet parcel s lisici se plochou: {}".format(count_dif))
# kolik parcel ma sikmou plochu rozdilnou od vodorovne o vice nez 0.01
with fiona.open("parcely.shp", 'r', driver="ESRI Shapefile",
crs="epsg:5514", schema=schema) as ds:
with rasterio.open('/tmp/LITC67.tif', 'r') as litc:
count = 0
count_dif = 0
for f in ds:
count += 1
try:
a = clip_raster(litc, shape(f['geometry']))
prumerne_prevyseni = slope_avg(a)
sikma_plocha = math.sqrt(
prumerne_prevyseni ** 2 + shape(f["geometry"]).area ** 2
)
rozdil = abs(sikma_plocha - shape(f["geometry"]).area)
if rozdil > tolerance:
count_dif += 1
except ValueError:
pass
print("-" * 80)
print("Počet parcel: {}".format(count))
print("Počet parcel s lisici se plochou: {}".format(count_dif))