Pollution Control by Identifying Potential Land for Afforestation - Python Project

Last Updated : 23 Jul, 2025

The program aims at controlling the pollution in a given area by suggesting the number of trees and the areas where they should be planted. The heart of the program is Computer Vision. A sample image is given below to get an idea about what we are going to do in this article. Note that we are going to implement this project using the Python language.

Tools and Technology used

In this project, we are using numpy and maths for calculation of our surrounding areas, PIL(Python Imaging Library) for manipulating. Before jumping to the project let’s understand these terms.

NumPy (Numerical Python): NumPy is a library for the Python programming language, adding support for large, multi-dimensional arrays and matrices, along with a large collection of high-level mathematical functions to operate on these arrays.
Maths: The Python math module offers you the ability to perform common and useful mathematical calculations within your application. Here are a few practical uses for the math module: Calculating combinations and permutations using factorials. Calculating the height of a pole using trigonometric functions.
PIL(Python Imaging Library): Python Imaging Library is a free and open-source additional library for the Python programming language that adds support for opening, manipulating, and saving many different image file formats.
OpenCV: OpenCV is a cross-platform library using which we can develop real-time computer vision applications. It mainly focuses on image processing, video capture and analysis including features like face detection and object detection.

Step by Step Implementation

Step 1: Create a New Project

Create a new project in PyCharm IDE or V.S. Code

Step 2: Before going to the coding section first you have to do some pre-task

In this project, we need an API key provided by Google Maps.

Step 3: Let's code Map segmentation

The satellite image generated by the 1st step undergoes Image segmentation, which separates all the objects in the image by focussing on edges and boundaries. The image is divided into objects such as buildings, trees, water bodies, roads, barren land, etc. Our first algorithm of choice is Mean Shift Algorithm for segmentation.

Python

import numpy as np
import cv2
from PIL import Image
import urllib.parse
import urllib.request
import io
from math import log, exp, tan, atan, pi, ceil
from place_lookup import find_coordinates
from calc_area import afforestation_area

EARTH_RADIUS = 6378137
EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS
INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0
ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0

def latlontopixels(lat, lon, zoom):
    mx = (lon * ORIGIN_SHIFT) / 180.0
    my = log(tan((90 + lat) * pi / 360.0)) / (pi / 180.0)
    my = (my * ORIGIN_SHIFT) / 180.0
    res = INITIAL_RESOLUTION / (2 ** zoom)
    px = (mx + ORIGIN_SHIFT) / res
    py = (my + ORIGIN_SHIFT) / res
    return px, py

def pixelstolatlon(px, py, zoom):
    res = INITIAL_RESOLUTION / (2 ** zoom)
    mx = px * res - ORIGIN_SHIFT
    my = py * res - ORIGIN_SHIFT
    lat = (my / ORIGIN_SHIFT) * 180.0
    lat = 180 / pi * (2 * atan(exp(lat * pi / 180.0)) - pi / 2.0)
    lon = (mx / ORIGIN_SHIFT) * 180.0
    return lat, lon

query = input('What kinda places you want me look up? ')
results = find_coordinates(query)

zoom = 18

ullat, ullon = results['upper_left']
lrlat, lrlon = results['lower_right']

scale = 1
maxsize = 640

ulx, uly = latlontopixels(ullat, ullon, zoom)
lrx, lry = latlontopixels(lrlat, lrlon, zoom)

dx, dy = lrx - ulx, uly - lry

cols, rows = int(ceil(dx / maxsize)), int(ceil(dy / maxsize))

bottom = 120
largura = int(ceil(dx / cols))
altura = int(ceil(dy / rows))
alturaplus = altura + bottom

final = Image.new("RGB", (int(dx), int(dy)))
for x in range(cols):
  
    for y in range(rows):
        dxn = largura * (0.5 + x)
        dyn = altura * (0.5 + y)
        latn, lonn = pixelstolatlon(ulx + dxn, uly - dyn - bottom / 2, zoom)
        position = ','.join((str(latn), str(lonn)))
        print(x, y, position)
        urlparams = urllib.parse.urlencode({'center': position,
                                            'zoom': str(zoom),
                                            'size': '%dx%d' % (largura, alturaplus),
                                            'maptype': 'satellite',
                                            'sensor': 'false',
                                            'scale': scale,
                                            'key': 'AIzaSyA_d4uV3HqPPWbCb77VhXNYn5UcXRLAiVc'})
        urlparamsmaps = urllib.parse.urlencode({'center': position,
                                                'zoom': str(zoom),
                                                'size': '%dx%d' % (largura, alturaplus),
                                                'maptype': 'roadmap',
                                                'sensor': 'false',
                                                'scale': scale,
                                                'key': 'AIzaSyA_d4uV3HqPPWbCb77VhXNYn5UcXRLAiVc'})
        url = 'http://maps.google.com/maps/api/staticmap?' + urlparams
        url1 = 'http://maps.google.com/maps/api/staticmap?' + urlparamsmaps
        f = urllib.request.urlopen(url)
        h = urllib.request.urlopen(url1)
        image = io.BytesIO(f.read())
        imagemaps = io.BytesIO(h.read())
        im = Image.open(image)
        immaps = Image.open(imagemaps)
        im.save("map.png")
        immaps.save("map_normal.png")

        img = cv2.imread('map.png')
        img_maps = cv2.imread('map_normal.png')
        shifted = cv2.pyrMeanShiftFiltering(img, 7, 30)
        shifted_normal = cv2.pyrMeanShiftFiltering(img_maps, 7, 30)
        gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
        ret, thresh = cv2.threshold(
            gray, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
        hsv = cv2.cvtColor(shifted, cv2.COLOR_BGR2HSV)
        hsv_normal = cv2.cvtColor(shifted_normal, cv2.COLOR_BGR2HSV)

        lower_trees = np.array([10, 0, 30])
        higher_trees = np.array([180, 100, 95])

        lower_houses = np.array([90, 10, 100])
        higher_houses = np.array([255, 255, 255])

        lower_roads = np.array([0, 0, 250])
        higher_roads = np.array([20, 20, 255])

        lower_feilds = np.array([0, 50, 100])
        higher_feilds = np.array([50, 255, 130])

        lower_feilds_blue = np.array([0, 80, 100])
        higher_feilds_blue = np.array([255, 250, 255])

        masktree = cv2.inRange(hsv, lower_trees, higher_trees)
        maskhouses = cv2.inRange(hsv, lower_houses, higher_houses)
        maskroads = cv2.inRange(hsv_normal, lower_roads, higher_roads)
        maskfeilds = cv2.inRange(hsv, lower_feilds, higher_feilds)
        gausssion_blur_maskfields = cv2.GaussianBlur(maskfeilds, (15, 15), 0)
        gausssion_blur_masktree = cv2.GaussianBlur(masktree, (15, 15), 0)
        blue_limiter = cv2.inRange(hsv, lower_feilds_blue, higher_feilds_blue)
        res_roads = cv2.bitwise_and(img_maps, img, mask=maskroads)
        # res_houses = cv2.bitwise_and(img,img,mask=maskhouses)
        res_feilds = cv2.bitwise_and(img, img, mask=gausssion_blur_maskfields)
        res_trees = cv2.bitwise_and(img, img, mask=masktree)

        # show the output image
        cv2.imshow('res', res_trees)
        cv2.imshow('res_fields', res_feilds)
        cv2.imshow('res_roads', res_roads)
        
        # cv2.imshow('res_houses',res_houses)
        # cv2.imshow('mask',maskfeilds)
        cv2.imshow('img', img)
        
        # cv2.imshow("hsv", hsv)
        cv2.waitKey(0)
        cv2.destroyAllWindows()

tot_land_area_acres, number_of_trees = afforestation_area()

Step 4: Let's code for place lookup

The user is required to input the name of the area, on which the program has to be executed. The satellite images of that area will be scraped and are zoomed such as to generate a clear image of the map on which image segmentation can be done.

Python

import urllib.parse
import requests

def find_coordinates(query):

    main_api = 'https://maps.googleapis.com/maps/api/place/textsearch/json?'
    url = main_api + \
        urllib.parse.urlencode({'query': query}) + '&key=Your API Key'

    json_data = requests.get(url).json()
    json_status = json_data['status']
    print('\nAPI Status :' + json_status)

    if json_status == 'OK':

        location_details = {
            'name_of_place': json_data['results'][0]['name'],
            'formatted_address': json_data['results'][0]['formatted_address'],
            'location': json_data['results'][0]['geometry']['location'],
            'upper_left': (json_data['results'][0]['geometry']['viewport']['northeast']['lat'],
                           json_data['results'][0]['geometry']['viewport']['southwest']['lng']),
            'lower_right': (json_data['results'][0]['geometry']['viewport']['southwest']['lat'], 
                            json_data['results'][0]['geometry']['viewport']['northeast']['lng']),
        }

        return location_details

Step 5: Let's code calculate area

We will be finding the pollution level of the given area. According to that level, we will find the number of trees required to bring that particular level to normal. In this process, we need to train a Classifier that can identify the buildings, the trees, and most importantly, the free land. The Zernike moments used by the above method will be used as features for these segments. The classifier is trained with labels as 'buildings', 'trees', 'water', 'free land', and 'roads'. After the training, we only need to find the part coming under the 'Free Land' label.

Python

import cv2
import numpy as np

def afforestation_area():

    img = cv2.imread('map.png')
    shifted = cv2.pyrMeanShiftFiltering(img, 7, 30)
    gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
    ret, thresh = cv2.threshold(
        gray, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
    hsv = cv2.cvtColor(shifted, cv2.COLOR_BGR2HSV)

    lower_trees = np.array([10, 0, 10])
    higher_trees = np.array([180, 180, 75])

    lower_houses = np.array([90, 10, 100])
    higher_houses = np.array([255, 255, 255])

    lower_roads = np.array([90, 10, 100])
    higher_roads = np.array([100, 100, 100])

    lower_feilds = np.array([0, 20, 100])
    higher_feilds = np.array([50, 255, 255])

    lower_feilds_blue = np.array([0, 80, 100])
    higher_feilds_blue = np.array([255, 250, 255])

    masktree = cv2.inRange(hsv, lower_trees, higher_trees)
    maskhouses = cv2.inRange(hsv, lower_houses, higher_houses)
    maskroads = cv2.inRange(hsv, lower_roads, higher_roads)
    maskfeilds_houses = cv2.inRange(hsv, lower_feilds, higher_feilds)
    blue_limiter = cv2.inRange(hsv, lower_feilds_blue, higher_feilds_blue)
    maskfeilds = maskfeilds_houses
    res = cv2.bitwise_and(img, img, mask=maskfeilds)

    print(res.shape)  # (640, 622, 3)
    print(np.count_nonzero(res))  # 679089

    print("number of pixels", res.size//3)
    tot_pixels = res.size//3
    # print("number of pixels: row x col", res.)

    no_of_non_zero_pixels_rgb = np.count_nonzero(res)
    row, col, channels = res.shape  # 152886
    print("percentage of free land : ", (no_of_non_zero_pixels_rgb /
                                         (row*col*channels)))  # 0.5686369573954984
    percentage_of_land = no_of_non_zero_pixels_rgb/(row*col*channels)

    # https://www.unitconverters.net/typography/centimeter-to-pixel-x.htm
    # says 1 cm = 37.795275591 pixels
    cm_2_pixel = 37.795275591
    print("row in cm ", row/cm_2_pixel)
    print("col in cm ", col/cm_2_pixel)

    row_cm = row/cm_2_pixel
    col_cm = col/cm_2_pixel
    tot_area_cm = tot_pixels/(row_cm*col_cm)
    tot_area_cm_land = tot_area_cm*percentage_of_land

    print("Total area in cm^2 : ", tot_area_cm_land)

    # in google maps 2.2cm = 50m => 1cm = 22.727272727272727m
    # in real life at zoom 18 1cm^2 = (22.727272727272727m)^2
    # = 516.5289256198347 m^2
    print("Total area in m^2 : ", tot_area_cm_land*(516.5289256198347))
    tot_area_m_actual_land = tot_area_cm_land*(516.5289256198347)

    # 1 m^2 = 0.000247105 acres :: source Google
    tot_area_acre_land = tot_area_m_actual_land*0.000247105
    print("Total area in acres : ", tot_area_acre_land)

    # https://www.treeplantation.com/tree-spacing-calculator.html
    # says if you have 2 ft between rows, and 2ft between 
    # trees will can take 10890 trees per acre.
    number_of_trees = tot_area_acre_land*10890
    print(f"{round(number_of_trees)} number of trees can be planted in\
    {tot_area_acre_land} acres.")

    return tot_area_acre_land, round(number_of_trees)

    # show the output image
    # cv2.imshow('res',res)

    # cv2.imshow('mask',maskfeilds)
    # cv2.imshow('img', img)

    #cv2.imshow("hsv", hsv)
    # cv2.waitKey(delay=0)
    # cv2.destroyAllWindows()

# afforestation_area()

Step 6: Let's code main file

Python

import numpy as np
import cv2
from PIL import Image
import urllib.parse
import urllib.request
import io
from math import log, exp, tan, atan, pi, ceil
from place_lookup import find_coordinates

EARTH_RADIUS = 6378137
EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS
INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0
ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0

def latlontopixels(lat, lon, zoom):
    mx = (lon * ORIGIN_SHIFT) / 180.0
    my = log(tan((90 + lat) * pi / 360.0)) / (pi / 180.0)
    my = (my * ORIGIN_SHIFT) / 180.0
    res = INITIAL_RESOLUTION / (2 ** zoom)
    px = (mx + ORIGIN_SHIFT) / res
    py = (my + ORIGIN_SHIFT) / res
    return px, py

def pixelstolatlon(px, py, zoom):
    res = INITIAL_RESOLUTION / (2 ** zoom)
    mx = px * res - ORIGIN_SHIFT
    my = py * res - ORIGIN_SHIFT
    lat = (my / ORIGIN_SHIFT) * 180.0
    lat = 180 / pi * (2 * atan(exp(lat * pi / 180.0)) - pi / 2.0)
    lon = (mx / ORIGIN_SHIFT) * 180.0
    return lat, lon

def calculate_area(res):
    """
    Args:
        Takes the transformed image as input
    Returns:
        :tot_area_acre_land: empty area in acres.
        :trees: rounded number of trees in the possible area.
    """
    # print(res.shape) # (640, 622, 3)
    # print(np.count_nonzero(res)) # 679089

    # print("number of pixels", res.size//3)
    tot_pixels = res.size//3
    # print("number of pixels: row x col", res.)

    no_of_non_zero_pixels_rgb = np.count_nonzero(res)
    row, col, channels = res.shape  # 152886
    
    percentage_of_land = no_of_non_zero_pixels_rgb/(row*col*channels)

    # https://www.unitconverters.net/typography/centimeter-to-pixel-x.htm
    # says 1 cm = 37.795275591 pixels
    cm_2_pixel = 37.795275591
    # print("row in cm ", row/cm_2_pixel)
    # print("col in cm ", col/cm_2_pixel)

    row_cm = row/cm_2_pixel
    col_cm = col/cm_2_pixel
    tot_area_cm = tot_pixels/(row_cm*col_cm)
    tot_area_cm_land = tot_area_cm*percentage_of_land

    # print("Total area in cm^2 : ", tot_area_cm_land)

    # in google maps 2.2cm = 50m => 1cm = 22.727272727272727 m 
    # in real life at zoom 18 1cm^2 = (22.727272727272727m)^2 
    # = 516.5289256198347 m^2
    tot_area_m_actual_land = tot_area_cm_land*(516.5289256198347)

    # 1 m^2 = 0.000247105 acres :: source Google
    tot_area_acre_land = tot_area_m_actual_land*0.000247105
    # print("Total area in acres : ", tot_area_acre_land)

    # https://www.treeplantation.com/tree-spacing-calculator.html
    # says if you have 2 ft between rows, and 2ft between trees 
    # will can take 10890 trees per acre.

    number_of_trees = tot_area_acre_land*10890
    # print(f"{round(number_of_trees)} number of trees can be planted
    # in {tot_area_acre_land} acres.")

    return tot_area_acre_land, round(number_of_trees)

def air_pollution_core(ullat, ullon, lrlat, lrlon, results):

    zoom = 18
    scale = 1
    maxsize = 640

    ulx, uly = latlontopixels(ullat, ullon, zoom)
    lrx, lry = latlontopixels(lrlat, lrlon, zoom)

    dx, dy = lrx - ulx, uly - lry

    cols, rows = int(ceil(dx / maxsize)), int(ceil(dy / maxsize))

    bottom = 120
    largura = int(ceil(dx / cols))
    altura = int(ceil(dy / rows))
    alturaplus = altura + bottom

    final = Image.new("RGB", (int(dx), int(dy)))
    total_acres_place, total_trees = 0., 0.
    total_tile_results = dict()
    for x in range(cols):
        for y in range(rows):
            dxn = largura * (0.5 + x)
            dyn = altura * (0.5 + y)
            latn, lonn = pixelstolatlon(
                ulx + dxn, uly - dyn - bottom / 2, zoom)
            position = ','.join((str(latn), str(lonn)))
            # print(x, y, position)
            urlparams = urllib.parse.urlencode({'center': position,
                                                'zoom': str(zoom),
                                                'size': '%dx%d' % (largura, alturaplus),
                                                'maptype': 'satellite',
                                                'sensor': 'false',
                                                'scale': scale,
                                                'key': 'YOUR_API_HERE'})
            url = 'http://maps.google.com/maps/api/staticmap?' + urlparams
            f = urllib.request.urlopen(url)
            image = io.BytesIO(f.read())
            im = Image.open(image)
            im.save("map_{}_{}_{}.png".format(x, y, position))

            img = cv2.imread("map_{}_{}_{}.png".format(x, y, position))
            shifted = cv2.pyrMeanShiftFiltering(img, 7, 30)
            gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
            ret, thresh = cv2.threshold(
                gray, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
            hsv = cv2.cvtColor(shifted, cv2.COLOR_BGR2HSV)

            lower_trees = np.array([10, 0, 10])
            higher_trees = np.array([180, 180, 75])

            lower_houses = np.array([90, 10, 100])
            higher_houses = np.array([255, 255, 255])

            lower_roads = np.array([90, 10, 100])
            higher_roads = np.array([100, 100, 100])

            lower_feilds = np.array([0, 20, 100])
            higher_feilds = np.array([50, 255, 255])

            lower_feilds_blue = np.array([0, 80, 100])
            higher_feilds_blue = np.array([255, 250, 255])

            masktree = cv2.inRange(hsv, lower_trees, higher_trees)
            maskhouses = cv2.inRange(hsv, lower_houses, higher_houses)
            maskroads = cv2.inRange(hsv, lower_roads, higher_roads)
            maskfeilds_houses = cv2.inRange(hsv, lower_feilds, higher_feilds)
            blue_limiter = cv2.inRange(
                hsv, lower_feilds_blue, higher_feilds_blue)
            maskfeilds = maskfeilds_houses
            res = cv2.bitwise_and(img, img, mask=maskfeilds)

            area_in_acres, number_of_trees = calculate_area(res)
            total_acres_place += area_in_acres
            total_trees += number_of_trees
            # print(f"area: {area_in_acres}, no of trees: {number_of_trees}")

            tile_results = {
                "name_of_tile_image": "map_{}_{}_{}.png".format(x, y, position),
                "area_acres": area_in_acres,
                "number_of_trees": number_of_trees
            }
            # print(tile_results)
            total_tile_results["{}_{}_{}".format(
                x, y, position)] = tile_results
            
            # uncomment below for viewing the output images
            # cv2.imshow('res',res)
            # cv2.imshow('img', img)
            # cv2.waitKey(delay=2000)
            # cv2.destroyAllWindows()
            
    # print(total_tile_results)
    results["total_tile_results"] = total_tile_results
    results["total_acres_of_land"] = total_acres_place
    results["total_number_of_trees"] = total_trees
    return results


def location_based_estimation(place):
    """
    :place: is a string that expects a name of a place
    """
    results = find_coordinates(place)

    ullat, ullon = results['upper_left']
    lrlat, lrlon = results['lower_right']

    returning_json = air_pollution_core(ullat, ullon, lrlat, lrlon, results)
    return returning_json


def coordinates_based_estimation(ullat, ullon, lrlat, lrlon):
    """
    :upperleft: a string expecting upperleft coordinates of 
    the tile you are expecting. ex : '12.92,79.11'
    :lowerright: a string expecting lowerright coordinates of
    the tile you are expecting. ex :'12.91,79.13'
    """
    # print(f"{upperleft.replace('\"','')}")
    # ullat, ullon = map(float, upperleft.split(','))
    # lrlat, lrlon = map(float, lowerright.split(','))
    results = dict()

    returning_json = air_pollution_core(ullat, ullon, lrlat, lrlon, results)
    return returning_json