2. crime location analysis by
   1. truncated negative exponential
   2. python3 source code, uses folium
3. 1. test program only, maybe buggy
5. here murder on Kyllikki Saari, Isojoki, Finland
7. python 3 script
8. 1. 30.8.2023 v. 0000.0000.0002

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 import scipy.stats as st

1. normal1=1/np.sqrt(2*np.pi)*np.exp(-np.power(x,2)/2)

dimx=500 dimy=500

kenterlon1=21.951671 kenterlat1=62.144996

1. olon1=21.92
2. olat1=62.125
3. olon2=21.98
4. olat2=62.155

olon1=21.88 olat1=62.08 olon2=22.00 olat2=62.20

1. param_f = 1
2. param_g = 2
3. dist_coeff=10

1. 1. for truncatednegative exponnetial distr

param_f = 1 ## inverse outer decay rate param_g = 2 ## inverse inner decay rate dist_coeff=1

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 ver 2 ...

xk=0 yk=0

for n in range (1, lenu): xk=lokx[n]+xk yk=loky[n]+yk

xk=xk/lenu yk=yk/lenu

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

xr=0 yr=0

for n in range (1, lenu): xr2=lokx[n]-xk yr2=loky[n]-yk xr=xr+xr2 yr=yr+yr2

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

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. ux=x-i
7. uy=y-j
8. if(ux==0): ux=1
9. if(uy==0): uy=1
10. distance=(math.sqrt(ux*ux+uy*uy))/buffer_dist
11. ## normal, ring radius 1,
12. p0=st.norm.pdf(distance,1,0.5)
13. #p0=math.log(math.pow(distance, -1))
15. p=p+p0
17. rosmo1[y,x]=p ## tsek it!
18. p = 0

for x in range(dimx): for y in range(dimy): p0=0 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 ## normal, ring radius 1, #p=st.norm.pdf(distance,param_f,param_g) p=0 #if(distance<0): # #dis2=1-distance # dis2=distance # p=math.exp(-dis2/param_g)

#p0=math.log(math.pow(distance, -1)) #if(p>p0): p0=p+p0 #p=math.exp(-distance) if(distance<=1): p=math.exp(-abs(1-distance)/param_g) if(distance>1): p=math.exp(-abs(1-distance)/param_f) p0=p0+p

rosmo1[y,x]=p0

plt.imshow(rosmo1)

plt.show()

1. quit(-1)

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 = ['black','black']

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=32

1. print()

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

1. geojson = geojsoncontour.contourf_to_geojson(contourf=contourf,min_angle_deg=0.01,ndigits=5,stroke_width=0.5,fill_opacity=0.5)

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)

1. Plot the contour on Folium map
2. folium.GeoJson(geojson,style_function=lambda x: {'color':"#f0f000f",'weight':0.5,'fillColor': None,'opacity': 0.5,}).add_to(map)

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

folium.raster_layers.ImageOverlay(rosmo2,

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

folium.LayerControl().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()