I try to make an A-star Path Finding program base on tile-graph (image raster data as graph) where each pixel values represent as cost and elevation.. In my case, DEM raster data was called by using gdal, then converting to numpy array to be processed on A* path finding implementation. This similar with Least Cost-Path analysis in GIS.
This is my standalone script i've been made so far :
import gdal
from gdalconst import *
import numpy as np
from numpy import savetxt
import math as mt
def Input_Dem():
filename = "F:\dem.tif"
dataset = gdal.Open( filename, GA_ReadOnly )
if dataset is None:
print "failed"
return dataset
def Input_Cost():
filename = "F:\cost.tif"
dataset = gdal.Open( filename, GA_ReadOnly )
if dataset is None:
print "failed"
return dataset
def Weight():
W = 1
return W
def RasterMatrix():
dem_surface = Input_Dem()
cost_surface = Input_Cost()
DEM_Matrix = np.array(dem_surface.GetRasterBand(1).ReadAsArray())
COST_Matrix = np.array(cost_surface.GetRasterBand(1).ReadAsArray())
geotransform_dem = dem_surface.GetGeoTransform()
geotransform_cost = cost_surface.GetGeoTransform()
col_dem = dem_surface.RasterXSize
row_dem = dem_surface.RasterYSize
col_cost = dem_surface.RasterXSize
row_cost = cost_surface.RasterYSize
cell_id_elevation = dict(enumerate(DEM_Matrix.flatten(),1))
cell_id_cost = dict(enumerate(COST_Matrix.flatten(),1))
cell_size = geotransform_dem[1]
col = col_dem
row = row_dem
return cell_id_elevation, cell_id_cost, col, row, cell_size, DEM_Matrix, COST_Matrix
def Var():
start_id = 1
end_id = 151044
open_list = [ ]
close_list = [ ]
current_cell_id = 0
accumulated_cost = { start_id : 0 }
parent_list = { }
heuristic_type = 'M'
return start_id, end_id, open_list, close_list, parent_list, accumulated_cost, current_cell_id, heuristic_type
def adjacent_cell(o,dx,dy):
adjacent_square = [o-(dx+1),o-dx,o-(dx-1),o-1,o+1,o+(dx-1),o+dx,o+(dx+1)]
for i in range(dx,dy*dx+dx,dx) :
if o == i :
adjacent_square = filter(lambda n: n != o-(dx-1) and n != o+1 and n != o+(dx+1) ,adjacent_square)
for i in range(1,dy*dx,dx) :
if o == i :
adjacent_square = filter(lambda n: n != o-(dx+1) and n != o-1 and n != o+(dx-1) ,adjacent_square)
adjacent_square = filter(lambda n: n > 0 and n <= dx*dy,adjacent_square)
return adjacent_square
def H(ht,y1,y2,x1,x2):
if ht == 'M' :
h = abs(x2 - x1) + abs(y2 - y1)
elif ht == 'E' :
h = mt.sqrt((abs(x2-x1))**2+(abs(y2-y1))**2)
elif ht == 'D' :
h = max(abs(x2-x1),abs(y2-y1))
return h
def Eucl_Dist(p,o,h,u):
P_Eucl_Dist = mt.sqrt(u**u + pow((h[p]-h[o]),2))
D_Eucl_Dist = mt.sqrt(2*u**u + pow((h[p]-h[o]),2))
return P_Eucl_Dist, D_Eucl_Dist
def Slope_W(p,o,h,u,w):
P_Slope_W = mt.degrees(mt.atan((h[p]-h[o])/u))
D_Slope_W = mt.degrees(mt.atan((h[p]-h[o])/mt.sqrt(2)*u))
return P_Slope_W, D_Slope_W
def F(p,o,h,u,w,c,dx,e_id,ht,CC):
y1 = mt.ceil(p/float(dx))
y2 = mt.ceil(e_id/float(dx))
x1 = p - (dx * (y1-1))
x2 = e_id - (dx * (y2-1))
if o - dx == p or o + dx == p or o + 1 == p or o - 1 == p :
G = ( (c[o] + c[p])/2. + Slope_W(p,o,h,u,w)[0] ) * Eucl_Dist(p,o,h,u)[0] + CC[o]
elif o-(dx+1) == p or o-(dx-1) == p or o+(dx-1) == p or o+(dx+1) == p :
G = ( (c[o] + c[p])/2. + Slope_W(p,o,h,u,w)[1] ) * Eucl_Dist(p,o,h,u)[1] + CC[o]
F = G + H(ht,y1,y2,x1,x2)
return F
def A_star(h,c,dx,dy,u,s_id,e_id,Op,Cl,Prt,CC,o,ht,w):
Op.append(s_id)
while e_id not in Op :
candidate = { }
for i in Op :
d = {i : CC[i]}
candidate.update(d)
o = min(candidate, key=candidate.get)
Cl.append(o)
Op.remove(o)
adjacent_list = adjacent_cell(o,dx,dy )
for p in adjacent_list :
if p in Cl:
adjacent_list = filter(lambda i: i != p, adjacent_list)
elif p not in Op : #
Op.append(p)
d = {p : o }
Prt.update(d)
d = {p : F(p,o,h,u,w,c,dx,e_id,ht,CC)}
CC.update(d)
elif id in Op :
f1 = F(p,o,h,u,w,c,dx,e_id,ht,CC)
f2 = F(p,Prt[p],h,u,w,c,dx,e_id,ht,CC)
if f1 < f2 :
d = {p : o }
Prt.update(d)
d = {id : F(p,o,h,u,w,c,dx,e_id,ht,CC)}
CC.update(d)
def root_path(Prt, e_id) :
yield e_id
while e_id in Prt:
e_id = Prt[e_id]
yield e_id
PATH_ID = list(root_path(Prt, e_id))
print PATH_ID
PATH_Matrix = [ ]
for i in xrange(1,dx*dy+1,1) :
if i in PATH_ID :
PATH_Matrix.append(1) # h[i]
else :
PATH_Matrix.append(0)
return np.array(PATH_Matrix).reshape(dy,dx)
path = A_star(RasterMatrix()[0],RasterMatrix()[1],RasterMatrix()[2],RasterMatrix()[3],RasterMatrix()[4],Var()[0],Var()[1],Var()[2],Var()[3],Var()[4],Var()[5],Var()[6],Var()[7],Weight())
print "path =", path
The result from this script is an array matrix representing a path connecting two location ( two cell ). For example (dim :25 x 29 ) :
path =
[[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]]
If i use a small resolution of DEM, 25x29 as an input, this program goes well. But if i try a large resolution of DEM, for example 1228 x 972. This program took very long time to calculate a whole cell. I guest that the problem is on the looping progress where i try to use loop for iteration (Line 110, A-star function).
while e_id not in Op :
Is there any solution to make my code run faster ?
