File:Log equi distance function of kyllikki saari murder site 1 v 1 r 1.png

2. 1. log distance function of points in mat
3. 1. test program only, maybe buggy
5. here murder on Kyllikki Saari, Isojoki, Finland
7. python 3 script
8. 1. 25.8.2023 v. 0000.0000.0000

import os import time import math from math import pow import numpy as np import matplotlib.pyplot as plt from PIL import Image

import folium from folium import plugins from selenium import webdriver

import branca from scipy.interpolate import griddata import geojsoncontour import scipy as sp import scipy.ndimage

dimx=512 dimy=512

kenterlon1=21.951671 kenterlat1=62.144996

olon1=21.92 olat1=62.125 olon2=21.98 olat2=62.155

param_f = 1 param_g = 2 dist_coeff=10

sw1=[ olat1, olon1] ne1=[ olat2, olon2]

rosmo1=np.zeros((dimy, dimx))

rosmo1=rosmo1*1

data =[[ 62.144996, 21.951671, 1. ],

      [    62.133412, 21.956548,   1.    ]]

1. [ 62.1406, 21.9481, 1. ]]

1. 1. minimum radius oh highest propab

1. 1. calculatend in below

1. buffer_dist = 0

data2=np.array(data)

dimlon=olon2-olon1 dimlat=olat2-olat1 lonk=dimx/dimlon latk=dimy/dimlat

lokx=[] loky=[]

print(olon1, olat1, olon2, olat2)

for n in range (0, len(data2)): lon0=data2[n,1] lat0=data2[n,0] print(lon0, lat0) x=int((lon0-olon1)*lonk) y=dimy-int((lat0-olat1)*latk) #x=int((lon0-olon1)*lonk) #y=int((lat0-olat1)*latk) lokx.append(x) loky.append(y)

print(lokx) print(loky)

lenu=len(data2)

sites=[]

for n in range (0, lenu): x=lokx[n] y=loky[n] #print(x,y) pp=(x,y) sites.append(pp)

print(sites)

1. 1. experimental only auto-radius ...

x00=lokx[0] y00=loky[0]

xr=0 yr=0

for n in range (1, lenu): x=lokx[n]-x00 y=loky[n]-y00 xr=xr+x yr=yr+y

xr=xr/lenu yr=yr/lenu

rr=int(math.sqrt(xr*xr+yr*yr))

print(rr)

buffer_dist =rr

1. 1. buffer_dist =0

1. 1. check this

1. quit(-1)

p = 0

def phi( x, y, i, j ):

   if abs(x-i) + abs(y-j) > buffer_dist:
       return 1
   else:
       return 0

def phi2( x, y, i, j ):

   if (math.sqrt( abs(x-i)*abs(x-i) + abs(y-j)*abs(y-j)) ) > buffer_dist:
       return 1
   else:
       return 0

1. 1. rossmo, manhattan distance

1. for x in range(dimx):
2. for y in range(dimy):
3. for site in sites:
4. i = site[0]
5. j = site[1]
6. characteristic = phi(x, y, i, j)
7. manhattan_dist = abs(x-i) + abs(y-j)
8. try:
9. p += characteristic / pow(manhattan_dist, param_f)
10. except ZeroDivisionError:
11. pass

1. n = (1-characteristic)*pow(buffer_dist, (param_g-param_f))
2. d = (2*buffer_dist)-characteristic
3. d = pow(d, param_g)
4. rosmo1[y,x]=p ## tsek it!
5. p = 0

   #print("")

1. 1. simple distance function

p=0

for x in range(dimx):

   for y in range(dimy):
       for site in sites:
           i = site[0]
           j = site[1]
           ux=x-i
           uy=y-j 
           if(ux==0): ux=1
           if(uy==0): uy=1       
           distance=(math.sqrt(ux*ux+uy*uy))/buffer_dist
           p0=math.log(math.pow(distance, -1))
           p=p+p0

       rosmo1[y,x]=p ## tsek it!
       p = 0

imr=np.copy(rosmo1) imb=np.copy(rosmo1) img=np.copy(rosmo1) ima=np.copy(rosmo1)

imr=imr*1 imb=imb*1 img=img*0 ima=255-imr

ima2=ima/256

rosmo2=np.dstack((imr, img, imb, ima))

1. rosmo2=np.power(rosmo2, 1.5)*10

deltalon=dimlon/dimx deltalat=dimlat/dimy

lontab=[] lattab=[] rostab=[]

for iy in range(0, dimy-1): tlat=olat2-deltalat*iy for ix in range(0, dimx-1): tlon=olon1+deltalon*ix tros=rosmo1[iy,ix] lontab.append(tlon) lattab.append(tlat) rostab.append(tros)

lontab=np.asarray(lontab) lattab=np.asarray(lattab) rostab=np.asarray(rostab)

print (np.shape(lontab)) print (np.shape(lattab)) print (np.shape(rostab))

1. quit(-1)

vmin=np.min(rostab) vmax=np.max(rostab)

lon_arr = np.linspace(np.min(lontab), np.max(lontab), dimx) lat_arr = np.linspace(np.min(lattab), np.max(lattab), dimy) lon_mesh, lat_mesh = np.meshgrid(lon_arr, lat_arr)

ros_mesh = griddata((lontab, lattab), rostab, (lon_mesh, lat_mesh), method='linear')

sigma = [5, 5]

1. ros_mesh = sp.ndimage.filters.gaussian_filter(ros_mesh, sigma, mode='constant')

levels=18

1. colors = ['blue','royalblue', 'navy','pink', 'mediumpurple', 'darkorchid', 'plum', 'm', 'mediumvioletred', 'palevioletred', 'crimson',
2. 'magenta','pink','red','yellow','orange', 'brown','green', 'darkgreen']

colors = ['blue','red']

levels = len(colors)

cm = branca.colormap.LinearColormap(colors, vmin=vmin, vmax=vmax).to_step(levels)

1. levels=10
2. levels=[1,10,100,200]

levels=1

1. print()

contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.75, colors=colors, linestyles='None', vmin=vmin, vmax=vmax)

geojson = geojsoncontour.contourf_to_geojson(

   contourf=contourf,
   min_angle_deg=0.01,
   ndigits=5,
   stroke_width=0.5,
   fill_opacity=0.1)

1. quit(-1)

map = folium.Map(location=[kenterlat1, kenterlon1], TileProvider='OpenStreetMap', zoom_start=16, control_scale = True,)

1. map = folium.Map(location=[kenterlat1, kenterlon1], TileProvider='Stamen Toner', zoom_start=16, control_scale = True,)
2. map = folium.Map(location=[kenterlat1, kenterlon1], TileProvider='stamentoner’', zoom_start=16)

1. plugins.HeatMap(data, radius = 50, min_opacity = 0.5, max_val = 1,gradient={.01: 'blue', 0.5: 'yellow', 1: 'red'}).add_to(map)

folium.raster_layers.ImageOverlay(rosmo2,

                   [[ olat1, olon1],[ olat2, olon2]],
                   opacity=0.5,
                  ).add_to(map)

folium.LayerControl().add_to(map)

1. Plot the contour on Folium map
2. folium.GeoJson(
3. geojson,
4. style_function=lambda x: {
5. 'color': "#ff0000",
6. 'weight': 0.1,
7. 'fillColor': None,
8. 'opacity': 0.15,
9. }
10. ).add_to(map)

1. folium.GeoJson( geojson, ).add_to(map)

1. 1. markers
   2. grave
2. folium.Marker(location=[data[0][0], data[0][1]],html="Hauta",popup="Kyllikki saaren suohauta").add_to(map)
   1. bike
3. folium.Marker(location=[data[1][0], data[1][1]],html="Pyörä", popup="Polkupyörän kätköpaikka").add_to(map)

latx1=data[0][0] lonx1=data[0][1]

folium.Circle([latx1, lonx1], 45, fill_color="grey", opacity=1.0,color = 'grey', fill=True).add_child(folium.Popup('Hauta')).add_to(map) folium.map.Marker(

   [latx1 + 0.0001,  lonx1 - 0.0001],
   icon=folium.DivIcon(
       icon_size=(150,36),
       icon_anchor=(0,0),

       )
   ).add_to(map)
   
  

latx1=data[1][0] lonx1=data[1][1]

folium.Circle([latx1, lonx1], 45, fill_color="blue", opacity=1.0,color = 'blue', fill=True).add_child(folium.Popup('Pyörä')).add_to(map) folium.map.Marker(

   [latx1 + 0.0001,  lonx1 - 0.0001],
   icon=folium.DivIcon(
       icon_size=(150,36),
       icon_anchor=(0,0),

       )
   ).add_to(map)
   

map.fit_bounds([sw1, ne1])

1. map.save("rossmotest1.html")

delay=5

1. Save the map as an HTML file

fn='testmap.html' tmpurl='file://{path}/{mapfile}'.format(path=os.getcwd(),mapfile=fn) map.save(fn)

1. quit(-1)

1. Open a browser window...

browser = webdriver.Firefox()

1. ..that displays the map...

browser.get(tmpurl)

1. Give the map tiles some time to load

time.sleep(delay)

1. Grab the screenshot

browser.save_screenshot('testmap.png')

1. Close the browser

browser.quit()