El programa tiene como objetivo controlar la contaminación en un área determinada al sugerir la cantidad de árboles y las áreas donde se deben plantar. El corazón del programa es Computer Vision. A continuación se muestra una imagen de muestra para tener una idea de lo que vamos a hacer en este artículo. Tenga en cuenta que vamos a implementar este proyecto utilizando el lenguaje Python .
Herramientas y Tecnología utilizadas
En este proyecto, estamos utilizando numpy y matemáticas para el cálculo de nuestras áreas circundantes, PIL (Biblioteca de imágenes de Python) para manipular. Antes de saltar al proyecto, entendamos estos términos.
- NumPy (Python numérico): NumPy es una biblioteca para el lenguaje de programación Python, que agrega soporte para arreglos y arrays grandes y multidimensionales, junto con una gran colección de funciones matemáticas de alto nivel para operar en estos arreglos.
- Matemáticas : el módulo matemático de Python le ofrece la posibilidad de realizar cálculos matemáticos comunes y útiles dentro de su aplicación. Aquí hay algunos usos prácticos para el módulo de matemáticas: Cálculo de combinaciones y permutaciones usando factoriales. Cálculo de la altura de un poste mediante funciones trigonométricas.
- PIL (Biblioteca de imágenes de Python): La Biblioteca de imágenes de Python es una biblioteca adicional gratuita y de código abierto para el lenguaje de programación Python que agrega soporte para abrir, manipular y guardar muchos formatos de archivo de imagen diferentes.
- OpenCV : OpenCV es una biblioteca multiplataforma con la que podemos desarrollar aplicaciones de visión artificial en tiempo real. Se enfoca principalmente en el procesamiento de imágenes, la captura y el análisis de videos, incluidas funciones como la detección de rostros y la detección de objetos.
Implementación paso a paso
Paso 1: Crear un nuevo proyecto
Cree un nuevo proyecto en PyCharm IDE o VS Code
Paso 2: antes de ir a la sección de codificación, primero debe hacer una tarea previa
En este proyecto, necesitamos una clave API proporcionada por Google Maps.
Paso 3: Codifiquemos la segmentación del mapa
La imagen de satélite generada por el primer paso se somete a la segmentación de la imagen, que separa todos los objetos de la imagen enfocándose en los bordes y los límites. La imagen se divide en objetos como edificios, árboles, masas de agua, carreteras, terrenos baldíos, etc. Nuestro primer algoritmo de elección es el algoritmo de desplazamiento medio para la segmentación.
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'})
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()
Paso 4: codifiquemos para la búsqueda de lugares
El usuario debe ingresar el nombre del área en la que debe ejecutarse el programa. Las imágenes satelitales de esa área se rasparán y se ampliarán para generar una imagen clara del mapa en el que se puede realizar la segmentación de imágenes.
Python3
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
Paso 5: Codifiquemos el área para calcular
Encontraremos el nivel de contaminación del área dada. De acuerdo con ese nivel, encontraremos la cantidad de árboles necesarios para llevar ese nivel particular a la normalidad. En este proceso, necesitamos entrenar a un clasificador que pueda identificar los edificios, los árboles y, lo que es más importante, el terreno libre. Los momentos de Zernike utilizados por el método anterior se utilizarán como características para estos segmentos. El clasificador está entrenado con etiquetas como ‘edificios’, ‘árboles’, ‘agua’, ‘tierra libre’ y ‘caminos’. Después de la capacitación, solo necesitamos encontrar la parte que se encuentra bajo la etiqueta ‘Tierra libre’.
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)
# 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()
Paso 6: Codifiquemos el archivo principal
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)
# 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
Producción:
Así es como se ve la estructura completa del proyecto.
Enlace GitHub: https://github.com/abhishektyagi2912/airpollution-analyses