-
Notifications
You must be signed in to change notification settings - Fork 0
/
BM2_AddRasterToVector.py
65 lines (53 loc) · 1.8 KB
/
BM2_AddRasterToVector.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
# Building Height Estimation
# Prepare input data ... Add Satellite imagery data to building footprint
# Thepchai Srinoi
# Department of Survey Engineering
import geopandas as gpd
import rasterio as rio
from rasterstats import zonal_stats
import pandas as pd
import numpy as np
from tqdm import tqdm
# https://gis.stackexchange.com/questions/433246/zonal-statistics-spatial-join-raster-to-polygon-in-python
# Building Footprint
zones = "building_inbound_DSM.gpkg"
# Satellite Image Raster data
values = "thesis_train_18082023.tif"
# I recommend .... DIVIDED AND CONQUER
# height consideration : more than this value
height_criteria = 20
#########################################################################
# Processing
# Import Geodataframe ... convert to UTM Zone 47
print('import geodataframe')
gdf = gpd.read_file(zones)
gdf = gdf.to_crs(32647)
#gdf = gdf.iloc[0:10]
print( gdf )
# Height Filtering
gdf = gdf.loc[ gdf.height > height_criteria]
gdf = gdf.reset_index(drop=True)
#print(gdf.head())
# Import DSM
print('import raster')
img = rio.open(values)
b_count = img.count
#import pdb; pdb.set_trace()
# Zonal Statistics ... Median Add Height to Vector Data
for myband in range(1,b_count+1) :
print('band ...' , myband )
print('rasterstat')
stats = gpd.GeoDataFrame(zonal_stats(gdf, values, stats=["median"]
, band= myband))
print('====================================')
print( stats )
print('add to dataframe')
#gdf['median'] = stats['median']
gdf[str(myband)] = stats
print( gdf )
#print( gdf.head() )
# Export Data
print('output')
gdf.to_file('traindata'+str(height_criteria) +'.gpkg', driver='GPKG')
gdf.to_csv('traindata'+str(height_criteria) +'.csv')
import pdb; pdb.set_trace()