Import libraries

In [1]:
import pandas as pd

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/**/*.tif'
DSM_tif = tif_file(path)
DSM_tif[:3]
Out[4]:
['./DSM/DHMVIIDSMRAS1m_k1/GeoTIFF/DHMVIIDSMRAS1m_k1.tif',
 './DSM/DHMVIIDSMRAS1m_k2/GeoTIFF/DHMVIIDSMRAS1m_k2.tif',
 './DSM/DHMVIIDSMRAS1m_k3/GeoTIFF/DHMVIIDSMRAS1m_k3.tif']
In [5]:
path = './DTM/**/*.tif'
DTM_tif = tif_file(path)
DTM_tif[:3]
Out[5]:
['./DTM/DHMVIIDTMRAS1m_k1/GeoTIFF/DHMVIIDTMRAS1m_k1.tif',
 './DTM/DHMVIIDTMRAS1m_k2/GeoTIFF/DHMVIIDTMRAS1m_k2.tif',
 './DTM/DHMVIIDTMRAS1m_k3/GeoTIFF/DHMVIIDTMRAS1m_k3.tif']

Geocoding - single address

In [6]:
address = 'Sint-Jorisstraat 93 8730 Beernem'
In [7]:
# enter an address in Belgium
# address = input("Enter an address in 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))
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Convert longtitude & latitude

In [8]:
# longtitude and latitude of the address entered
lat,lon = lat_long(address)
lat,lon
Out[8]:
(51.1273793, 3.3386131)
In [9]:
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 [10]:
# coordinate EPSG:31370
xx,yy = transform(lon, lat)

longtitude:77886.85955147157 latitude:202423.4265439231

Using CSV to locate coordinate EPSG:31370

In [11]:
path = '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses_Wallonie/Belgium_Address.csv'
In [12]:
def get_coordinate(path,address):
    
    index = []
    
    df = pd.read_csv(path)
    
    for i in range(df.shape[0]):

        if df['Address'][i] == address :
            index.append(i)
            break
        continue

    xx, yy = df['EPSG:31370_x'][index][index[0]] ,df['EPSG:31370_y'][index][index[0]] 

    lat, lon = df['EPSG:4326_lat'][index][index[0]] ,df['EPSG:4326_lon'][index][index[0]] 
    
    return (xx, yy) , (lat, lon)
In [13]:
(xx, yy) , (lat, lon) = get_coordinate(path,address)

Locate tif for given longtitude & latitude

In [14]:
# 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 [15]:
# 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/DHMVIIDSMRAS1m_k'+ str(i) +'/GeoTIFF/DHMVIIDSMRAS1m_k'+ str(i) + '.tif')  

            print('The house is in this tif :', 'DHMVIIDSMRAS1m_k' + str(i) + '.tif')
            
    return scr_path

dsm_path = check_DSM(xx,yy)

The house is in this tif : DHMVIIDSMRAS1m_k13.tif
In [16]:
# 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/DHMVIIDTMRAS1m_k'+ str(i) +'/GeoTIFF/DHMVIIDTMRAS1m_k'+ str(i) + '.tif')  

            print('The house is in this tif :', 'DHMVIIDTMRAS1m_k' + str(i) + '.tif')
            
    return scr_path

dtm_path = check_DTM(xx,yy)

The house is in this tif : DHMVIIDTMRAS1m_k13.tif

Clip the location of the house from tif

In [17]:
dsm_path
Out[17]:
['./DSM/DHMVIIDSMRAS1m_k13/GeoTIFF/DHMVIIDSMRAS1m_k13.tif']
In [18]:
dtm_path
Out[18]:
['./DTM/DHMVIIDTMRAS1m_k13/GeoTIFF/DHMVIIDTMRAS1m_k13.tif']
In [23]:
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(address+"_clipped_dsm.tif",tiled=True, dtype="int32")

    return clipped.plot()
In [24]:
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(address +"_clipped_dtm.tif",tiled=True, dtype="int32")

    return clipped.plot()
In [25]:
dsm_path[0]
Out[25]:
'./DSM/DHMVIIDSMRAS1m_k13/GeoTIFF/DHMVIIDSMRAS1m_k13.tif'
In [26]:
clip_dsm(dsm_path[0],20)
Out[26]:
<matplotlib.collections.QuadMesh at 0x7f947587eb50>
In [27]:
clip_dtm(dtm_path[0],20)
Out[27]:
<matplotlib.collections.QuadMesh at 0x7f948cfb4850>
In [28]:
def chm_tif():
    
    with rasterio.open(address +'_clipped_dtm.tif') as src:
        lidar_dtm_im = src.read(1, masked=True)
        dtm_meta = src.profile
    
    with rasterio.open(address +'_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(address +'_clipped_chm.tif', 'w', **dsm_meta) as ff:
        ff.write(lidar_chm,1)
        
    chm_tif = address +'_clipped_chm.tif'
    
    return chm_tif

Create 3D

In [29]:
def House_3D(tif,size:int,azim=215):

    chm = imageio.imread(tif)

    nx,ny = chm.shape

    x = np.linspace(0, size*2, nx)
    y = np.linspace(0, size*2, ny)

    yv,xv = 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=azim)
        
    fig.colorbar(chm3d, shrink=0.3, aspect=10);
    fig.savefig(address +'_3D.png', dpi=200) 
    
    return  plt.show()
In [30]:
House_3D(chm_tif(),20,45)
In [32]:
# compare google map image of the house

from IPython.display import Image
Image(address+".png")
Out[32]:
In [31]:
# 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[31]:
True
In [33]:
## Plotting the CHM with Mayavi2
 
from osgeo import gdal
from mayavi import mlab

ds = gdal.Open(chm_tif())
data = ds.ReadAsArray()

data = data.astype(np.float32)

mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

surf = mlab.surf(data) 
mlab.show()
In [2]:
Image(address+"_mayavi.png")
Out[2]: