Model for green area detection in cities — TUTORIAL for the City of Buenos Aires

In previous posts, we analyzed the importance of green area detection to promote quality public policies and take care of the environment. Moreover, we reviewed the main steps for mapping them, through the use of machine learning tools and satellite imagery: pre-processing, model training, prediction and post-processing

In this post, we explain each phase in detail and we show you step by step how to detect green areas in cities and classify them applying a classification model based on convolutional neural networks.

At the end of this tutorial, we will have generated a map of green areas in the City of Buenos Aires (our example case), classified in two classes: parks and vegetation (public woodland, boulevards and wastelands).

Green areas detected and classified as parks (left) or vegetation (right).

Downloading the data

We will start by downloading the datasets of interest. In this tutorial, we will use the open-source tools Unetseg and Satproc, aerial images of the City of Buenos Aires and annotations of parks and vegetation in the city.

• Downloading aerial images of Buenos Aires neighborhoods (it includes Belgrano, Chacarita, Colegiales, Palermo, Recoleta, Villa Crespo).

import os

from glob import glob

barrio='Belgrano'

path_to_images= './workshop/imagenes/'

!mkdir -p $path_to_images # Creates the directory

bucket_imagen_caba=os.path.join('gs://dym-workshops-public/gcba/workshop/imagen_2021/2021_RGB_utm.tif') # Image of CABA used for training

!gsutil -m cp -r $bucket_imagen_caba $path_to_images # Downloads the image

• Downloading the vector file containing annotations: this file contains polygons of the classes parks and vegetation.

path_to_gt=’./workshop/data/gt/’
!mkdir -p $path_to_gt
bucket_gt=’gs://dym-workshops-public/gcba/workshop/anotaciones/areas_verdes.gpkg’
!gsutil -m cp -r $bucket_gt $path_to_gt

• Downloading the vector file of the area of interest: the vector file of the area of interest defines the areas where there are annotations. By defining them, we reduce the computing time, since only part of the image is analyzed.

path_to_aoi=’./workshop/data/shp/ ‘
!mkdir -p $path_to_aoi
bucket_aoi=’gs://dym-workshops-public/gcba/workshop/data/shp/aoi_areas_verdes.geojson’
!gsutil -m cp -r $bucket_aoi $path_to_aoi

Pre-processing: Satproc tool

The purpose of the pre-processing step is to generate the training dataset based on an adequate number of high-quality images in the required format. The segmentation models require a dataset of images and masks which delimit the object of interest. These are binary images: 1 indicates where there is an object, and 0 indicates where there are no objects.

To achieve that, we will use our tool Satproc, which precisely allows us to create a dataset of images and masks in the format required by some machine learning models of image segmentation.

How does it work? Satproc analyses each image through a search window of a fixed size (size). If it finds a polygon within the image, it generates a mask with pixel values, where 1 indicates that there is a polygon, and 0 indicates that there are no polygons. Then, it moves forward to the next step (step_size) and it repeats the process. If more than one class is identified, it concatenates the masks in a single window, generating an n-bands image, where n indicates the number of classes.

Next, it takes the image and the vector file with the annotations and generates a directory with two folders: one containing the images dataset and another one containing the masks dataset (one mask per image). Since our goal is to detect two classes of objects, the mask will be a binary image of two bands: one for vegetation and one for parks.

Example image (left) and its masks: vegetation mask (middle) and park mask (right).

1.1 Generating the training dataset

As mentioned before, to generate the training dataset we need satellite imagery with the band combination we want to use, and a vector file with annotations of the ground truth polygons (.geojson .shp .gpkg, etc), both with the same georeference.

For this analysis, we used EPSG:32721. Besides, the vector file must include a “class” field identifying each polygon class, as it is shown in the figure.

Class
1 P
2 V
3 V
4 V
5 P
6 V

> Install Satproc and some other libraries that we will use:

!pip install pysatproc

> To generate the images for the training dataset:

size = 3000
step_size = 1000
path_dataset_train=os.path.join(“./workshop/dataset/data_train/” +str(size)+’_’+str(step_size)+’/’)
!mkdir -p $path_dataset_trainvector_file = ‘./workshop/data/gt/areas_verdes.gpkg’
aoi= “./workshop/data/shp//aoi_areas_verdes.geojson”
path_to_files = ‘./workshop/imagenes/2021_RGB_utm.tif’ # Image of CABA
output_folder = os.path.join(“./workshop/dataset/data_train/” +str(size)+’_’+str(step_size)+’/’)
!satproc_extract_chips \
$path_to_files \
-o $output_folder \
 — size $size \
 — step-size $step_size \
 — aoi $aoi \
 — labels $vector_file \
 — label-property ‘class’ \
 — classes ‘V’ ‘P’

These are the arguments:

• o: destination path > we recommend a descriptive path, for example, “data_train/3000_1000/” describes the following: Data_train → data used for the training; 3000_1000 → (the images are square images of 3000x3000 and the step_size is 1000).
• size: size of the resulting images (the images are square).
• step-size: step of the process. If the step-size is the same as the size, there is no overlap in the images.
• Sometimes, generating chips with an overlap is useful for the training, since it provides us with more data to train. However, during the prediction, the value must be the same as the image size.
• aoi: path of the vector file where the neighborhoods are defined. By defining an area of interest, only the images from those neighborhoods are processed.
• Labels: vector file containing annotations.
• label-property: name of the field where each category is defined (it is only used for training).
• classes: names of the classes (as they appear in the .gpkg), divided by spaces.

This command will generate two folders in the destination path: “images” and “extent”. The type of the files of the images folder will be 3-bands Tiff (rgb) and the type of the files of the extent folder will also be Tiff, but with N bands, where N represents the number of classes. In this case, there are 2 bands, each of which is a binary mask.

2. Model training

This model is a deep neural network for segmenting CNN objects with a U-Net architecture. Originally, it was developed to segment biomedical imaging, but it has been found that it performs well in many other areas.

The architecture includes many down-sampling and up-sampling layers, which are connected. In this way, the information from the layers of greater accuracy activations is combined with the nodes of the previous layers of higher resolution.

U-net architecture

• Input: masks and images dataset.
• Output: for each image, the model generates an image where the value of each pixel represents the probability of finding the object of interest. If the prediction is made based on several classes, the output image will be an n-bands image, where each band represents each class.

The advantage of this kind of analysis is that the images may not be just RGB “visuals”. Other bands and georeferenced information, such as GIS layers, can be included to improve the results.

> First, we install Unetseg

This package streamlines the training and prediction of deep learning models based on the U-Net architecture for semantic segmentation problems in relation with geospatial imagery.

!pip install folium==0.2.1
!pip install unetseg
from unetseg.train import TrainConfig, train
from unetseg.evaluate import plot_data_generator
import os

> Then, we set the training:

Obs: It is useful to establish a name including information about the training parameters for the weights file (model name). For example: < images dim > < number of bands > < size >_< stepsize > < step_per_epoch > .h5 or similar.

model_name=’UNet_256x256_3N_2C_spe92_areas_verdes_modelo1.h5'

> Let’s determine how many images and masks includes the training dataset to calculate the parameter value: steps_per_epoch.

from glob import glob
path_to_images=glob(os.path.join(‘../workshop/dataset/data_train/’ +str(size)+’_’+str(step_size)+’/’+’images/*.tif’))
#path_to_mask=glob(os.path.join(“../workshop/dataset/data_train/” +str(size)+’_’+str(step_size)+’/’+’extent/*.tif’))
n=(len(path_to_images),len(path_to_mask))
print(n)
config = TrainConfig(width=256, # Size of the image processed by the UNet (it needs to be a multiple of 16, for example, 160, 320, etc; and higher than 80)
height=256,
n_channels=3, # Number of image channels, rgb -> 3
n_classes=2, # Number of classes which need to be classified
apply_image_augmentation=True, # If it is True, the dataset is expanded, generating new images based on small variations of the existing images (rotation),
seed=42,
epochs=20, # Number of times that the entire dataset can undergo the training process
batch_size=16, # Number of data that are processed each time, which can be limited by the gpu memory available (it must be a multiple of 16)
steps_per_epoch=92, # Generally, it has to be the same as the number of images / the batch_size; if it is higher, the number of images generated with image augmentation will increase
early_stopping_patience=5, # Throughout the training process, the results are saved after each epoch. If the error does not change after ¿¿10 ?? iterations, the process is interrupted because it is assumed that the error was already reduced significantly
validation_split=0.2, # The sample is split during training (80% of the data) and validation (20% of the data) to calculate the error during the training process
test_split=0.1, # Number of images of the dataset
images_path=os.path.join(“../workshop/dataset/data_train/” +str(size)+’_’+str(step_size)+’/’),# Path to the images generated with Satproc
 model_path=os.path.join(‘../data/weights/’, model_name),# Path of the model trained
model_architecture=’unet’,
evaluate=True,
class_weights= [0.2, 0.8])

Plot_data_generator allows the visualization of some images and masks of the training dataset.

plot_data_generator(num_samples=3, fig_size=(10, 10), train_config=config,img_ch = 2)

> Now, we run the training

res_config = train(config)

Metrics of the model

One of the most used metrics in computer vision is Intersection over Union (IoU).

Through this metric, the area of intersection between the ground truth and the prediction about the common area between both is calculated. The value tends to be 1 as the prediction improves.

Estimation of the Intersection over Union (IoU).

Estimation of the Intersection over Union (IoU)

import matplotlib.pyplot as plt
plt.figure(figsize=(16,4))
plt.subplot(121)
plt.plot(res_config.history[‘loss’])
plt.plot(res_config.history[‘val_loss’])
plt.title(‘Loss’)
plt.ylabel(‘Loss’)
plt.xlabel(‘Epoch’)
plt.legend([‘Train’, ‘Val’], loc=’upper left’)
plt.subplot(122)
plt.plot(res_config.history[‘mean_iou’])
plt.plot(res_config.history[‘val_mean_iou’])
plt.title(‘mean_iou’)
plt.ylabel(‘val_mean_iou’)
plt.xlabel(‘Epoch’)
plt.legend([‘Train’, ‘Val’], loc=’upper left’)
plt.show()

3) Prediction

Now, we will create the prediction dataset. Based on it, the model we trained before will generate predictions about these new images, determining the probability of finding objects of interest. To streamline this process, we will use Unetseg.

!pip install folium==0.2.1
!pip uninstall -y h5py
!pip install h5py==2.10.0
!pip install unetseg
!pip install pysatproc
from unetseg.predict import PredictConfig, predict
from unetseg.evaluate import plot_data_results
import os

> First, we need to download the image of the neighborhood

barrio=’Belgrano’
path_to_images=os.path.join(‘./workshop/images_predict/’+barrio+’/’)
!mkdir -p $path_to_images
bucket_imagen_barrio=os.path.join(‘gs://dym-workshops-public/gcba/workshop/’+barrio+’/imagen/’)
!gsutil -m cp -r $bucket_imagen_barrio $path_to_images

> Then, we will download the vector file corresponding to the area of interest

bucket_aoi_barrio=os.path.join(‘gs://dym-workshops-public/gcba/workshop/’+barrio+’/shp/’+barrio+’_1.geojson’)
aoi_predict=’./workshop/aoi_predict/’
!mkdir $aoi_predict
!gsutil -m cp -r $bucket_aoi_barrio $aoi_predict
> Finally, we will download the already trained model
path_to_model=’./workshop/data/’
!mkdir -p path_to_model
bucket_modelo=’gs://dym-workshops-public/gcba/workshop/modelo’
!gsutil -m cp -r $bucket_modelo $path_to_model

Now, we are ready to generate the prediction dataset. Similarly, we use satproc_extract_chips to generate the prediction dataset. In this case, the vector file defines the areas of interest for the prediction, and the size and the step_size are the same, since an image overlap is no longer useful.

size = 3000
step_size=size
path_to_files = os.path.join(‘./workshop/images_predict/’+barrio+’/imagen/*.tif’) # Image of the neighborhood
output_folder = os.path.join(‘./workshop/’+barrio+’/dataset/data_predict/’+str(size)+’_’+str(step_size)+’/’)
vector_file_aoi = os.path.join(‘./workshop/aoi_predict/’+barrio+’_1.geojson’)
!satproc_extract_chips \
$path_to_files \
-o $output_folder \
 — size $size \
 — step-size $step_size \
 — aoi $vector_file_aoi

> We check the name of the model

!ls ./workshop/data/

Next, we determine the prediction settings. And we need to transfer the path to the images of the prediction dataset and the model.

predict_config = PredictConfig(
images_path=os.path.join(‘./workshop/’+barrio+’/dataset/data_predict/’+str(size)+’_’+str(step_size)+’/’), # Path to the images on which the prediction will be based
results_path=os.path.join(‘./workshop/’+barrio+’/dataset/data_results/’+str(size)+’_’+str(step_size)+’/’), # Destination path for our prediction
batch_size=16,
model_path=os.path.join(‘./workshop/data/’,’UNet_256x256_3N_2C_spe92_areas_verdes_modelo1.h5'), # Path to the model (.h5)
height=256,
width=256,
n_channels=3,
n_classes=2,
class_weights= [0.2,0.8])
predict(predict_config) # Runs the prediction
plot_data_results(num_samples=2, fig_size=(5, 5), predict_config=predict_config, img_ch =2,n_bands =3)

4) Post-processing of the results

Finally, we have reached the post-processing step, which aims to improve the quality of the results. During this step, a threshold is applied to maintain the predictions with higher probabilities. In turn, once the final result is polygonized, a minimum area filter can be applied to reduce the number of “noisy” predictions.

The following feature applies a polygonization routine to the results of the model prediction and generates a vector file in GeoPackage format. (GPKG). The routine applies gdal_polygonize.py to each resulting chip, generating a GPKG for each of them, and then it efficiently combines all of these files into a single file.

Before combining them, it also applies a threshold to the raster values, which, in this case, represent the probability (values between 0 and 1). filter_by_max_prob is a Satproc feature that shows just the results which do not include any pixel exceeding the specified threshold.

> Now, we specify a probability threshold for each class:

threshold_A=50
threshold_B=50

The following code applies a probability threshold to each class and assigns the following values:

• A value of 100 to the pixels corresponding to the “Park” class.
• A value of 200 to the pixels corresponding to the “Vegetation” class.
• A value of 1 to the background.

The aim of this process is to convert the 2-bands probability image (binary) into a single-band image with those pixel values.

from tqdm import tqdm
import subprocess
import os
DEBUG=True
results_path = os.path.join(‘./workshop/’+barrio+’/dataset/data_results/’+str(size)+’_’+str(step_size)+’/’) # Path to the prediction results
# If DEBUG: print (results_path)
dest = os.path.join(‘./workshop/’+barrio+’/dataset/data_results/temp_postprocesamiento’,str(size)+”_”+str(step_size)) # Destination path for the post-processing
os.makedirs(dest, exist_ok=True)
#if DEBUG: print (dest)
files = next(os.walk(results_path))[2] # Analyses the files
new_dir = os.path.join(dest, ‘1band_2classes’,’thr_’+str(threshold_A)+’_’+str(threshold_B)) # We create a directory for the single-band image
os.makedirs(new_dir, exist_ok=True) # If the directory does not exist, one is created
#if DEBUG: print (new_dir)
for file in files:
if os.path.splitext(file)[-1]==’xml’: continue
result_dst = os.path.join(results_path,file)
output_tif_path = os.path.join(dest, ‘1band_2classes/’,’thr_’+str(threshold_A)+’_’+str(threshold_B),file) # Destination of the single-band image
exp = “((B<”+str(threshold_B)+”)*(A >”+str(threshold_A)+”)*199 ) + ((A<”+str(threshold_A)+”)*(B >”+str(threshold_B)+”)*99 ) +1" # Estimation to apply thresholds and assign values to the pixels, depending on their class
cmd_calc = f’gdal_calc.py -A {result_dst} — A_band=1 -B {result_dst} — B_band=2 — outfile {output_tif_path} — calc=”{exp}” — NoDataValue=0'
subprocess.run(cmd_calc, shell=True)
print(“ — “,output_tif_path)

> We convert the 2-bands probability image into a single-band image:

output_vrt = os.path.join(dest, ‘1band_2classes’,’thr_’+str(threshold_A)+’_’+str(threshold_B),’tiles.vrt’)
output_tif = os.path.join(dest, ‘1band_2classes’,’thr_’+str(threshold_A)+’_’+str(threshold_B),’*.tif’)
cmd_buildvrt = f’gdalbuildvrt {output_vrt} {output_tif}’
subprocess.run(cmd_buildvrt, shell=True)

> We execute the polygonization:

output_shp_path = os.path.join(‘./workshp/’+barrio+’/dataset/data_results/’, ‘poligonize_classAB_thr_’+str(threshold_A)+’_’+str(threshold_B)+’_’+str(size)+’_’+str(step_size)+’.gpkg’)
cmd_polygonize = f’gdal_polygonize.py {output_vrt} {output_shp_path}’
subprocess.run(cmd_polygonize, shell=True)

Optionally, we may apply a minimum area filter, but first we need to convert the results to meters (utm):

#src_file = ‘file_4326.gpkg’
#dst_file = ‘file_utm_32721.gpkg’
#!ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:32721 -f ‘GPKG’ $dst_file $src_file

> Finally, we filter the polygons which do not have the specified minimum area value:

min_area = 10
input_path = dst_file
output_path = “file_utm_32721.gpkg”
!ogr2ogr \
-t_srs EPSG:32721 \
-f “GPKG” \
-sql “SELECT * FROM tabla m WHERE (ST_Area(geom) > $min_area) “ \
-dialect SQLITE \
-nln results \
$output_path \
$input_path

If you follow these steps, you should obtain images like these:

Detection and classification of green areas: vegetation (red) and parks (blue)

You can access Google Colab to run the code from there and obtain more information about the detection of green areas in cities.

If you have any questions or if you would like to hear more about our solutions, follow us on our social networks or contact us.