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.
Python3
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'})
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.
Python3
import urllib.parse
import requests
def find_coordinates(query):
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.
Python3
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)
# 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)
# 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
Python3
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)
# 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)
# 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'})
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
|
Output:
This is how the complete project structure looks like.
GitHub Link: https://github.com/abhishektyagi2912/airpollution-analyses