Import libraries¶
In [1]:
from glob import glob
from natsort import natsorted, ns
import geopy
import folium
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter
import rasterio
from rasterio.plot import show
from pyproj import Transformer
from rasterio.windows import Window
import rioxarray
import imageio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
Search all tifs and sort files¶
In [2]:
def tif_file(path):
tif_file =[]
# using glob library to get all the file with .tif
files = glob(path,recursive = True)
for file in files:
tif_file.append(file)
# sort files with number in the file (ascending order) using natsorted library
tif_file = natsorted(tif_file, alg=ns.IGNORECASE)
return tif_file
In [3]:
# pwd '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses'
In [4]:
path = './DSM_PROV_BRABANT_WALLON/*.tif'
DSM_tif = tif_file(path)
DSM_tif
Out[4]:
In [5]:
path = './DTM_PROV_BRABANT_WALLON/*.tif'
DTM_tif = tif_file(path)
DTM_tif
Out[5]:
Geocoding - single address¶
In [6]:
# enter an address in Belgium
address = input("Enter an address in PROV_BRABANT_WALLON , Belgium:")
def lat_long(address):
# to get the longtitude and latitude of the address entered
locator = Nominatim(user_agent="myGeocoder")
location = locator.geocode(address)
house_lat_long = [location.latitude, location.longitude]
return house_lat_long
def house_map(func):
# to plot address on a map
house_map = folium.Map(location=func,zoom_start=150)
folium.Marker(location=func,popup=lat_long(address)).add_to(house_map)
house_map.save("house_map.html")
return house_map
house_map(lat_long(address))
# Example address:
# Rue Demaret 36 1300 Wavre
Out[6]:
Convert longtitude & latitude¶
In [7]:
# longtitude and latitude of the address entered
lat,lon = lat_long(address)
lat,lon
Out[7]:
In [8]:
def transform(lon, lat):
# transform to Belgium 'EPSG:31370' coordinate
transformer = Transformer.from_crs("EPSG:4326", crs_to = 'EPSG:31370' ,always_xy=True)
xx, yy = transformer.transform(lon, lat)
print(f'longtitude:{xx} latitude:{yy}' )
return xx,yy
In [9]:
# coordinate EPSG:31370
xx,yy = transform(lon, lat)
Locate tif for given longtitude & latitude¶
In [10]:
# create all bounding box from tifs
def bounds(tifs):
bounds = []
for i in tifs:
src = rasterio.open(i)
bounds.append(src.bounds)
return bounds
src_bounds = bounds(DSM_tif)
In [11]:
src_bounds[:3]
Out[11]:
In [12]:
# check which tif contains the location in DSM
def check_DSM(xx,yy):
scr_path = []
for i,bound in enumerate(src_bounds,1):
if (xx >= bound[0] and xx <= bound[2]) & \
(yy >= bound[1] and yy <= bound[3]):
scr_path.append('./DSM_PROV_BRABANT_WALLON/RELIEF_WALLONIE_MNS_2013_2014.tif')
print('The house is in this tif :', 'RELIEF_WALLONIE_MNS_2013_2014.tif')
else:
print('The house is NOT in this tif')
return scr_path
dsm_path = check_DSM(xx,yy)
In [13]:
# check which tif contains the location in DTM
def check_DTM(xx,yy):
scr_path = []
for i,bound in enumerate(src_bounds,1):
if (xx >= bound[0] and xx <= bound[2]) & \
(yy >= bound[1] and yy <= bound[3]):
scr_path.append('./DTM_PROV_BRABANT_WALLON/RELIEF_WALLONIE_MNT_2013_2014.tif')
print('The house is in this tif :', 'RELIEF_WALLONIE_MNT_2013_2014.tif')
else:
print('The house is NOT in this tif')
return scr_path
dtm_path = check_DTM(xx,yy)
Clip the location of the house from tif¶
In [14]:
dsm_path
Out[14]:
In [15]:
dtm_path
Out[15]:
In [16]:
def clip_dsm(path,window_size:int):
xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
# set window size
n = window_size
# create coordinates and geometries
coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]
coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]
geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]
# clip the image as per the geometries size
clipped = xds.rio.clip(geometries)
clip_dsm = clipped.rio.to_raster("clipped_dsm.tif",tiled=True, dtype="int32")
return clipped.plot()
In [17]:
def clip_dtm(path,window_size:int):
xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
# set window size
n = window_size
# create coordinates and geometries
coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]
coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]
geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]
# clip the image as per the geometries size
clipped = xds.rio.clip(geometries)
clip_dtm = clipped.rio.to_raster("clipped_dtm.tif",tiled=True, dtype="int32")
return clipped.plot()
In [18]:
dsm_path[0]
Out[18]:
In [19]:
clip_dsm(dsm_path[0],20)
Out[19]:
In [20]:
clip_dtm(dtm_path[0],20)
Out[20]:
In [21]:
def chm_tif():
with rasterio.open('clipped_dtm.tif') as src:
lidar_dtm_im = src.read(1, masked=True)
dtm_meta = src.profile
with rasterio.open('clipped_dsm.tif') as src:
lidar_dsm_im = src.read(1, masked=True)
dsm_meta = src.profile
lidar_chm = lidar_dsm_im - lidar_dtm_im
with rasterio.open('clipped_chm.tif', 'w', **dsm_meta) as ff:
ff.write(lidar_chm,1)
chm_tif = 'clipped_chm.tif'
return chm_tif
Create 3D¶
In [32]:
def House_3D(tif,size:int):
chm = imageio.imread(tif)
ny, nx = chm.shape
x = np.linspace(0, size*2, nx)
y = np.linspace(0, size*2, ny)
xv, yv = np.meshgrid(x, y)
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, projection='3d')
chm3d=ax.plot_surface(xv,yv,chm,cmap='plasma',linewidth=0)
ax.set_title('CHM(Canopy Height Model)')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Distance (m)')
ax.set_zlabel('Elevation (m)')
ax.view_init(azim=215)
fig.colorbar(chm3d, shrink=0.3, aspect=10);
fig.savefig('house.png', dpi=200)
return plt.show()
In [33]:
House_3D(chm_tif(),20)
In [24]:
# compare google map image of the house
from IPython.display import Image
Image("Rue Demaret 36.png")
Out[24]:
In [25]:
# to view actual house on google map
import webbrowser
url = 'https://www.google.com.my/maps/place/'+str(lat)+','+str(lon)
webbrowser.open(url)
Out[25]: