dynamics.py

from shapely.geometry import Polygon, Point, MultiPolygon
import folium
import pyproj
import json, os, random
import pandas as pd
from datetime import datetime

'''
This code aims to create a heatmap using records for the livestock. So the field is segmented with
a grid where the user can interact. 

Thus is as a ultimate script, because it use all data generated by the others script, i.e.,
'tcf_scraping.py' which builds the dataset about the vegetation, 'auravant_api.py' provides the 
field's NDVI, 'mongodb.py' feeds of BASTO dataset through MongoDB and returns the records of livestock
position and the fields limits.
'''

class HeatMap(object):
    '''<plots> is going to allocate the limits of the selected farm, and <heads> the historical
        record about the position for each cow with GPS'''
    plots, heads = 0, 0
    gps_path, plot_path = 0, 0

    ''' The object receives the farm's name and a float to control the grip density'''

    def __init__(self, farm: str, density: float):
        self.square_size = density

        self.farm = farm.replace(' ', '_')
        self.gps_path = './basto_dataset/gps_' + self.farm + '/'
        self.plot_path = './basto_dataset/plots/'
        try:
            ''' These lines are commented because the project deployed cannot reach the Mongo
                Database. However to run it locally it is suggested uncomment it. '''
            files = os.listdir(self.gps_path)
            self.heads = [pd.read_csv(self.gps_path + name) for name in files]
            for df in self.heads:
                df['timestamp'] = pd.to_datetime(df['timestamp'])

            with open(self.plot_path + self.farm + '.json', 'r') as f:
                self.plots = json.load(f)
        except:
            print('Get dataset by running: \n \
                  \t python3 mongodb.py')
            
    ''' From Auravant the limits of a field are brought as a string, e.g.,
            POLYGON(lng1 lat1, lng2 lat2, ... ...)
        Because of that, the following function transforms this string into the format
        the library shapely can deal with. For instance
         = [(lat1, lng1) , (lat2, lng2), ...]'''
    
    def vertex(self, polygon: str):
        this_vertex = []
        polygon = polygon.strip("POLYGON()").split(',')
        for p in polygon:
            p = p.split(' ')
            v = (float(p[1]), float(p[0]))
            this_vertex.append(v)
        return this_vertex
    
    '''<minimum_regions> scans all fields provided by the user, and take the union of all the fields
        so if there're fields into another field these are sweep it, and the algorith the most
        relevant regions to taking into account'''
    
    def minimum_regions(self):
        polygons = []
        for key in self.plots.keys():
            polygons.append(self.vertex_format(self.plots[key]))
        
        polygons = [Polygon(vertex) for vertex in polygons]
        union = polygons[0]
        for pp in polygons[1:]:
            try:
                union = union.union(pp)
            except:
                continue
        
        '''Finally it returns than regions which does not intersect to each other'''

        minimum = []
        for geom in union.geoms:
            minimum.append(list(geom.exterior.coords))
        
        return minimum
    
    def vertex_format(self, coords: dict):
        '''This transforms list of coordinates into shapely required format'''
        this_vertex = [(lat, lng) for lat, lng in zip(coords['lat'], coords['lng'])]
        return this_vertex
    
    '''<grid> method creates the grid on the provided <vertex> polygon'''
    
    def grid(self, polygon: list):
        '''This is the shapaly polygon object'''
        field = Polygon(polygon)
        '''field.envelope returns the coordinates of the minimun square that  contains the
            previous polygon'''
        square = field.envelope

        minx, miny, maxx, maxy = square.bounds

        ''' First, the previous square is divides into smaller squares according to the paramater
            <density>'''
        small_squares = []
        x = minx
        while x < maxx:
            y = miny
            while y < maxy:
                small_square_coords = [(x, y),
                                    (x + self.square_size, y),
                                    (x + self.square_size, y + self.square_size),
                                    (x, y + self.square_size),
                                    (x, y)]
                small_square = Polygon(small_square_coords)
                if small_square.intersects(square):
                    small_squares.append(small_square_coords)
                y += self.square_size
            x += self.square_size

        '''Finally it is chosen the regions intersecting the field'''
        intersection = []
        for sq in small_squares:
            solid = Polygon(sq)
            if solid.intersects(field):
                try:
                    inter = solid.intersection(field)
                    intersection.append(list(inter.exterior.coords))
                except:
                    for geom in inter.geoms:
                        intersection.append(list(geom.exterior.coords))

        return intersection
    
    '''<repoint> method uses the records located in heads and returns the multiple coordinates all
        GPS recorded over a time-range'''
    
    def repoint(self, from_date: datetime, to_date: datetime):
        position = []
        for points in self.heads:
            this_point = []
            df = points.loc[(points['timestamp'].dt.date >= pd.to_datetime(from_date).date()) & \
                            (points['timestamp'].dt.date <= pd.to_datetime(to_date).date())]
            for lat, lng in zip(df['lat'].values, df['lng'].values):
                this_point.append(Point(lat, lng))
            position.append(this_point)
        '''This returns a list of lists. Each element contains Point object of shapely library which
        depicts the coordinates of at certain time'''
        position = [x for x in position if x]
        return position
    
    '''<cmap> method builds two range, first to color out each cell on fucntions the number of points
        there's in the cell, namely, gives to each cell a number from 0 to 1, where 0 depicts
        the animal never was there, and 1 it always was there. And second to determinate the percet of
        time the animal spend time in the cell'''
    
    def cmap(self, intersection: list, points: list):
        this_cmap, this_cmap2 = {}, {}
        for ii, inter in enumerate(intersection):
            inter = Polygon(inter)
            this_cmap[ii] = 0
            for dot in points:
                if inter.contains(dot):
                    this_cmap[ii] += 1

        suma = sum([x for x in this_cmap.values()])

        for k, v in this_cmap.items():
            this_cmap2[k] = round(v / suma * 100, 1)

        cmin = min(this_cmap.values())
        cmax = max(this_cmap.values())

        for k, v in this_cmap.items():
            this_cmap[k] = round((v - cmin) / (cmax - cmin), 2)

        return this_cmap, this_cmap2
    
    '''The following method is the principal function, i.e., where the whole map is created.'''
    
    def heat_map(self, from_date: datetime, to_date: datetime):
        vertices = self.minimum_regions()
        intersection = []
        for minn in vertices:
            inter = self.grid(minn)
            for ii in inter:
                intersection.append(ii)

        polygon = [Polygon(vertex) for vertex in vertices]
        multipolygon = MultiPolygon(polygon)

        '''This points out the place where the map will be displayed'''
        centroid = multipolygon.centroid
        lat_center = centroid.coords[0][0]
        lon_center = centroid.coords[0][1]
        
        map = folium.Map(location=[lat_center, lon_center], zoom_start=13)

        '''As it can be appreciated there's two <points> parameters. The first one has the real data
            of the position of each animal in the livestock, it is commended to display the second one.
            This is formed by testing (non-real data), which is created randomly.
            
            So if you want to visualize the real date uncomment the first one and comment the socond one.
            You could comment <num_days> and <delta> as well.'''

        # points = self.repoint(from_date, to_date)
   

        delta = to_date - from_date
        num_days = delta.days
        points = self.test_points(multipolygon, num_days)
        
        for vertex in vertices:
            folium.Polygon(locations=vertex, color='green', fill=False, opacity=1).add_to(map)

        for kk, dots in enumerate(points):
            cmap, cmap2 = self.cmap(intersection, dots)
            self.paint(intersection, cmap, cmap2, map, kk)

        return map
    
    '''The following function receives a zone with its small regions, the cmaps and the map object,
        there's an id paramaeter but this is just a counter, however it can be replaces by the
        real GPS id. 
        With that paramater, the functions create the zones with more animals and so on'''
    
    def paint(self, intersection: list, cmap: dict, cmap2: dict, m: folium.Map, id: int):
        utm_proj = pyproj.Proj(proj='utm', zone=20, south=False)

        for ii, domain in enumerate(intersection):
            if cmap[ii] > 0.0:
                coords_m = [utm_proj(long, lat) for lat, long in domain]
                coords_m = Polygon(coords_m)
                area = round(coords_m.area / 10000, 3)

                polygon = folium.Polygon(locations=domain, color='white', fill=True, fill_color='red', opacity=0.4, fill_opacity=cmap[ii])
                popup = folium.Popup(f' _id: {id} <br> Area: {area} ha <br> Time: {cmap2[ii]}%', max_width=400)
                polygon.add_child(popup)
                m.add_child(polygon)

    '''This funcion has no relation to the main project but create a non-real dataset to display'''
    def test_points(self, multipolygon: MultiPolygon, days: int):
        minx, miny, maxx, maxy = multipolygon.bounds
        cenx, ceny = 0.001, 0.001
        varx, vary = days * cenx, days * ceny
        X, Y = [], []
        n = 10
        while len(X) != n:
            x = random.uniform(minx, maxx)
            y = random.uniform(miny, maxy)
            point = Point(x, y)
            if multipolygon.contains(point):
                X.append(x)
                Y.append(y)

        points = []
        m = n * 6
        for x, y in zip(X, Y):
            new_points = []
            while len(new_points) != m:
                x_ = x + random.normalvariate(cenx, varx)
                y_ = y + random.normalvariate(ceny, vary)
                point = Point(x_, y_)
                if multipolygon.contains(point):
                    new_points.append(point)
            points.append(new_points)
        
        return points