import random from shapely.geometry import Point, LineString, Polygon import matplotlib.pyplot as plt import geopandas as gpd import shapely import numpy as np # from scipy.spatial import Delaunay # from itertools import permutations from math import hypot ## No need to worry about these dummy = [] def randomPoints(polygon, min_dist, max_dist, pointsCount=3): assert min_dist < max_dist # make sure parameters are valid global points points = [] filled = False counter = 0 minx = polygon.bounds['minx'].min() miny = polygon.bounds['miny'].min() maxx = polygon.bounds['maxx'].max() maxy = polygon.bounds['maxy'].max() while len(points) < pointsCount: pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy)) if pnt.intersects(polygon.unary_union): pt_xy = (pnt.x, pnt.y) counter += 1 if counter > 10000: raise ValueError("Distance settings are too restrictive. Redefine the distance variables") if not filled: points.append(pnt) pointsAll.append(pt_xy) dummy.append(pnt) filled = True else: check_valid = 0 for item in points: if min_dist < item.distance(pnt) < max_dist: check_valid += 1 else: break if check_valid == len(points): points.append(pnt) pointsAll.append(pt_xy) dummy.append(pnt) return points def distanceCal(points_list, x, y): return points_list[x].distance(points_list[y]) def grids(ip_shape, grid_size): ip_poly = gpd.GeoDataFrame(geometry=[ip_shape], crs='epsg:3857', index=[0]) minx = ip_poly.bounds['minx'].min() miny = ip_poly.bounds['miny'].min() maxx = ip_poly.bounds['maxx'].max() maxy = ip_poly.bounds['maxy'].max() length = wide = grid_size cols = list(np.arange(minx, maxx + wide, wide)) rows = list(np.arange(miny, maxy + length, length)) polygons = [] for x in cols[:-1]: for y in rows[:-1]: polygons.append( Polygon([(x, y), (x+wide, y), (x+wide, y+length), (x, y+length)])) grid = gpd.GeoDataFrame({'geometry': polygons}, crs='epsg:3857') grid_clip = gpd.clip(grid, ip_poly) grid_clip['UID'] = range(len(grid_clip)) # Don't use *.index grid_clip['Area'] = grid_clip.area/ 10**6 area_mean = grid_clip["Area"].mean() grid_clip['flag'] = "" for row in grid_clip['Area'].index: if grid_clip.loc[row, "Area"] < (0.1*area_mean): grid_clip.loc[row, "flag"] = "Small area" return grid_clip ## Providing input ip_poly = 'POLYGON((7511352.944174381 7112117.215453198,7511653.217271762 7108850.245242323,7515833.018061557 7109162.529626478,7515556.765888274 7112537.598119422,7511352.944174381 7112117.215453198))' ip_shape = shapely.wkt.loads(ip_poly) polygonAll = grids(ip_shape, grid_size) finalize = False len(polygonAll) ## Making the points for SHP pointsAll = [] for i in range(len(polygonAll)): polygon = polygonAll.loc[polygonAll.UID == i] loc_area = polygon["Area"].values[0] ## Min. distance = 1/4 of the grid length local_min = np.sqrt(loc_area * 10**6)/4 randomPoints(polygon, min_dist= local_min, max_dist= 50000) ## Distance-matrix data = {} data['locations'] = pointsAll def compute_euclidean_distance_matrix(locations): """Creates callback to return distance between points.""" distances = {} for from_counter, from_node in enumerate(locations): distances[from_counter] = {} for to_counter, to_node in enumerate(locations): if from_counter == to_counter: distances[from_counter][to_counter] = 0 else: # Euclidean distance distances[from_counter][to_counter] = (int( hypot((from_node[0] - to_node[0]), (from_node[1] - to_node[1])))) return distances distance_matrix = compute_euclidean_distance_matrix(data['locations']) ## Adding the dummy point at the end with zero distance to every point dict_dummy_pt = {} for num in range(len(pointsAll)+1): dict_dummy_pt[num] = 0 distance_matrix[len(pointsAll)] = dict_dummy_pt for num in range(len(pointsAll)+1): distance_matrix[num][len(pointsAll)] = 0 ## o-r Tools from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp # import math def create_data_model(): """Stores the data for the problem.""" data['num_vehicles'] = 1 data['depot'] = 0 return data def print_solution(manager, routing, solution): """Prints solution on console.""" print('Objective: {} miles'.format(solution.ObjectiveValue())) index = routing.Start(0) plan_output = 'Route for vehicle 0:\n' route_distance = 0 while not routing.IsEnd(index): plan_output += ' {} ->'.format(manager.IndexToNode(index)) previous_index = index index = solution.Value(routing.NextVar(index)) route_distance += routing.GetArcCostForVehicle(previous_index, index, 0) plan_output += ' {}\n'.format(manager.IndexToNode(index)) print(plan_output) plan_output += 'Route distance: {}miles\n'.format(route_distance) def main(): """Entry point of the program.""" # Instantiate the data problem. data = create_data_model() distance_matrix = compute_euclidean_distance_matrix(data['locations']) # Create the routing index manager. manager = pywrapcp.RoutingIndexManager(len(distance_matrix), data['num_vehicles'], data['depot']) # Create Routing Model. routing = pywrapcp.RoutingModel(manager) def distance_callback(from_index, to_index): """Returns the distance between the two nodes.""" # Convert from routing variable Index to distance matrix NodeIndex. from_node = manager.IndexToNode(from_index) to_node = manager.IndexToNode(to_index) return distance_matrix[from_node][to_node] transit_callback_index = routing.RegisterTransitCallback(distance_callback) # Define cost of each arc. routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) # Setting first solution heuristic. search_parameters = pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy = ( routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC) # Solve the problem. solution = routing.SolveWithParameters(search_parameters) # # Print solution on console. # if solution: # print_solution(manager, routing, solution) def get_routes(solution, routing, manager): """Get vehicle routes from a solution and store them in an array.""" # Get vehicle routes and store them in a two dimensional array whose # i,j entry is the jth location visited by vehicle i along its route. routes = [] for route_nbr in range(routing.vehicles()): index = routing.Start(route_nbr) route = [manager.IndexToNode(index)] while not routing.IsEnd(index): index = solution.Value(routing.NextVar(index)) route.append(manager.IndexToNode(index)) routes.append(route) return routes routes = get_routes (solution, routing, manager) # Display the routes. # for i, route in enumerate(routes): # print('Route', i, route) # x =0 # # print(routes) global my_route my_route = routes[0][:-1] # print(my_route) if __name__ == '__main__': main() ## Making the route plottable as a map best_route_pt_list = [] for loc_pt in my_route: best_route_pt_list.append(dummy[int(loc_pt)]) best_route_pt_list from shapely.geometry import LineString route_polyline = LineString([z for z in best_route_pt_list]) route_geo = gpd.GeoSeries(route_polyline) print(best_route_pt_list[0]) ## Plotting as a map import matplotlib.pyplot as plt xs = [point[0] for point in pointsAll] ys = [point[1] for point in pointsAll] ax = polygonAll.plot(figsize=(20,20)) poly_boundary = polygonAll.boundary poly_boundary.plot(ax=ax, color='darkgray') print(route_geo.length) route_geo.plot(ax=ax, color='yellow') plt.scatter(xs, ys, c='black', s=50)