2. 1. 16.9.2023 v. 0000.0000.0021
4. 1. attempt to make geographical rossmo analysis of crime locations
   2. tries to shrikk area to find suspect of crime series
7. atlanta 1979 or jack the ripper
9. based on code at
10. 1. # https://medium.com/@errazkim/how-to-catch-serial-killers-using-one-simple-formula-29b1cbd58963
14. python 3 script

import os import time import math from math import pow import numpy as np

import scipy as sp import scipy.ndimage import scipy.stats as st from scipy.interpolate import griddata from scipy.spatial.distance import cdist from scipy.stats import norm from scipy.interpolate.rbf import Rbf

import matplotlib.pyplot as plt from matplotlib import ticker from matplotlib.ticker import StrMethodFormatter

from PIL import Image import netCDF4

import folium from folium import plugins from selenium import webdriver

import branca

import geojsoncontour

bufferi=0 ## auto set f=0.500 g=0.667

1. 1. settings !

1. rosmo matrix dims

dimx=100 dimy=100

1. 1. map initial zoom factor

zoomfac1=16

1. 1. draw

draw_suspect=1

1. 1. relative drawing skala of symbols

draw_skala=1

1. 1. extent of map

deltarectlon=0.12 deltarectlat=0.12

1. 1. symbols drawing scala

drawskala=1

1. 1. draw rossmo or another matrix

drawmatrix=1

1. 1. draw user image

drawimage=0

1. 1. actually is image
   2. https://commons.wikimedia.org/wiki/File:Jack_the_ripper_rossmo_3.png#/media/File:Whitechapel_Spitalfields_7_murders.JPG

imagename1="ripper1.png"

1. 1. corners of image https://commons.wikimedia.org/wiki/File:Jack_the_ripper_rossmo_3.png#/media/File:Whitechapel_Spitalfields_7_murders.JPG

1. 1. crime locations, Atlanta murders 1979-1981

crimelocations = np.array([[33.703093, -84.532406], [33.660032, -84.49509], [33.741141, -84.383959], [33.711061, -84.447227], [33.701493, -84.584169], [33.746652, -84.350482], [33.6605585, -84.4941276], [33.755227, -84.465294], [33.692281, -84.350066], [33.738875, -84.408613], [33.805397, -84.470401], [33.677783, -84.427292], [33.679033, -84.358048], [33.858629, -84.455166], [33.68205, -84.573247], [33.68164, -84.067554], [33.837342, -84.332364], [33.680747, -84.249387], [33.631259, -84.128966], [33.683782, -84.64159], [33.653495, -84.681008], [33.766465, -84.422132], [33.653852, -84.6803], [33.802901, -84.500141], [33.7392226, -84.3287373], [33.804113, -84.499154]])

1. 1. suspect locations, available in some cases
   2. Atlanta Williams

suspectlat = 33.754065 suspectlon = -84.446577 suspect = [suspectlat, suspectlon]

1. 1. jack the ripper

1. 1. ripper = np.array([[33.703093, -84.532406], [33.660032, -84.49509], [33.741141, -84.383959], [33.711061, -84.447227], [33.701493, -84.584169], [33.746652, -84.350482], [33.6605585, -84.4941276], [33.755227, -84.465294], [33.692281, -84.350066], [33.738875, -84.408613], [33.805397, -84.470401], [33.677783, -84.427292], [33.679033, -84.358048], [33.858629, -84.455166], [33.68205, -84.573247], [33.68164, -84.067554], [33.837342, -84.332364], [33.680747, -84.249387], [33.631259, -84.128966], [33.683782, -84.64159], [33.653495, -84.681008], [33.766465, -84.422132], [33.653852, -84.6803], [33.802901, -84.500141], [33.7392226, -84.3287373], [33.804113, -84.499154]])

1. 1. williams = [33.754065, -84.446577]

1. 1. crimenumbers=[1,2,3,4,5]
2. lats=[51.5203785, 51.5186270, 51.5199207, 51.5137500,51.5137633]
3. lons=[-0.0737000,-0.0752304,-0.0610239,-0.0779247,-0.0657004]
4. names=["Nichols","Eddowes","Kelly","Chapman","Stride"]
5. dates=['1888-01-01', '1888-01-01','1888-01-01','1888-01-01','1888-01-01']
6. ages=[43,46,35,47,44]

1. 1. jack the ripper
   2. nw
2. imglonnw1=-0.08173
3. imglatnw1=51.52137
4. imglonsw1=-0.08202
5. imglatsw1=51.51208

1. imglonne1=-0.05690
2. imglatne1=51.52137

1. imglonse1=-0.05679705085
2. imglatse1=51.51195544111
3. imglonsw1=-0.08202
4. imglatsw1=51.51208

1. imglonne1=-0.05690
2. imglatne1=51.52137

1. crimelocations=[[51.5203785,-0.0737000],[51.5186270,-0.0752304],[ 51.5199207, -0.0610239],[ 51.5137500,-0.0779247],[ 51.5137633, -0.0657004]]

1. suspectlat = 51.5203785
2. suspectlon = -0.0737000

1. suspect= [suspectlat, suspectlon]

suspectx = 0 suspecty = 0

1. 1. some next globals will be auto set bu script ...

numpoints=0

rectlon1=0 rectlon2=0 rectlat1=0 rectlat2=0 rectlatc=0 rectlatc=0 latsize=0 lonsize=0 rectmedianlon=0 rectmedianlat=0

dlat=0 dlon=0

kenterlat=0 kenterlon=0 kenterx=0 kentery=0

meanradius=0 devradius=0 meanpointsdistance=0 meanpointsdistance_manhattan=0 minimumdistance=0

medianx=0 mediany=0 medianlat=0 medianlon=0

minlat=0 maxlat=0 minlon=0 maxlon=0 cutlat=0 cutlon=0

suspectx = 0 suspecty = 0

koordis000=np.copy(crimelocations) koordis=np.copy(crimelocations)

lons=[] lats=[]

data=[]

datax=[] datay=[]

def crime_locations_dropper(selekt1): ## uses lon, lat coords global crimelocations, koordis global numpoints

crimelocations0=np.copy(crimelocations) numpoints0=len(crimelocations0)-1 #print(numpoints0)

cenx=0 ceny=0 rx=0 ry=0

for n in range(0, numpoints0): rax=crimelocations0[n][0] ray=crimelocations0[n][1] rx=rx+rax ry=ry+ray

cenx=rx/numpoints0 ceny=ry/numpoints0

rtab0=[] utab0=[]

for n in range(0, numpoints0): x=crimelocations0[n][0] y=crimelocations0[n][1] rax=crimelocations0[n][0]-cenx ray=crimelocations0[n][1]-ceny rr=math.sqrt(rax*rax+ray*ray) rtab0.append(rr) utab0.append(np.array([rr,x,y]))

#print(utab0)

utab1=np.array(utab0)

#print (np.shape(utab0))

utab2=utab1[utab1[:, 0].argsort()]

#print (np.shape(utab1)) #print(utab2)

#selekt1=6

utab3=utab2[:selekt1, :]

#print(utab3) utab4=utab3[:,1:] #print(utab4)

crimelocations=None

crimelocations=utab4 koordis=np.copy(crimelocations)

return(0)

1. 1. third rossmo

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

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

def rossmo_3(sites, budi, param_f, param_g): global dimx, dimy mata1=np.zeros((dimy, dimx)) p=0

for x in range(dimx): for y in range(dimy): p=0 for site in sites: i = site[0] j = site[1] characteristic = phi(budi, x, y, i, j) manhattan_dist = abs(x-i) + abs(y-j) try: p += characteristic / pow(manhattan_dist, param_f) except ZeroDivisionError: pass n = (1-characteristic)*pow(budi, (param_g-param_f)) d = (2*budi)-characteristic d = pow(d, param_g)

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

#mata1=np.fliplr(mata1) #mata1=np.fliplr(mata1) #mata1=np.rot90(mata1) mata1=np.flipud(mata1) return(mata1)

1. 1. second rossmo
2. def rossmo_formula_2(k: int, f: int, g: int, b: int, crimes:List[Tuple[int]], imax: int, jmax: int) -> Dict[str, Any]:

def rossmo_2(k,f,g, b, crimess): ## b buffer

global dimx, dimy imax=dimx jmax=dimy matra1 = np.zeros((imax, jmax), dtype=int) step = 1 j = 0 while (j < jmax): i=0 while (i < imax): crimes_sum = 0 for crime in (crimess): distance = abs(j - crime[0]) + abs(i - crime[1])

if distance > b: crimes_sum += 1/math.pow(distance, f) else: crimes_sum += math.pow(b, g-f)/math.pow(2*b - distance, g)

pij = k * crimes_sum matra1[i][j] = pij i += step j += step

matra1=np.flipud(matra1) return (matra1)

1. 1. first rossmo

def manhattan(A, B):

   return abs(A[0] - B[0]) + abs(A[1] - B[1])

def nearest(A, L):

   n = len(L)
   m = A
   d = 0
   p = 0
   j = 0
   while True:
     if A[0] != L[j][0] or A[1] != L[j][1]:
       #m = L[j]
       d = manhattan(A,L[j])
       p = j
       break 
     else: j += 1
   for i in range(len(L)):
     dis = manhattan(A,L[i])
     if dis != 0 and dis < d:
       d = dis
       p = i
   return L[p]

def radius_buffer(L):

   s = 0
   n = len(L)
   for i in range(n):
     M = L.copy()
     del M[i]
     B = nearest(L[i], M)
     s += manhattan(L[i],B)
   return s/(2*n)

def proba(i, j, L):

   # gives proba of criminal living in (i,j)
   # L = list of crime sites   (n,2) list
   proba = 0
   f = 1/3
   g = 2/3
   B = radius_buffer(L)
   for p in range(len(L)):
     d = manhattan(L[p], [i,j])
     if d > B:
       proba += 1/(d**f)
     else:
       proba += B**(g-f)/((2*B - d)**g)
   return proba

def rossmo_matrix(m,n,L): print("Generating rossmo matrix ...") M = [[0 for j in range(n)] for i in range(m)] for i in range(m): for j in range(n): M[i][j] += proba(i,j,L)

mata1=np.fliplr(M) mata1=np.rot90(mata1) mata1=np.flipud(mata1) return (mata1)

def generate_coords(): ## warning here negative koordis global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius global medianx,mediany, medianlat,medianlon global rectmedianlon, rectmedianlat, devradius global minlat, maxlat, minlon, maxlon, cutlon, cutlat

lons=koordis[:,1] lats=koordis[:,0]

numpoints=len(lons)

#print(lons) minlon1=np.min(lons) maxlon1=np.max(lons) minlat1=np.min(lats) maxlat1=np.max(lats) #print(minlon1)

rectmedianlon=(minlon1+maxlon1)/2 rectmedianlat=(minlat1+maxlat1)/2

kenterlon1=0 kenterlat1=0

for n in range(0, numpoints): kenterlat1=kenterlat1+lats[n] kenterlon1=kenterlon1+lons[n]

kenterlon=kenterlon1/numpoints kenterlat=kenterlat1/numpoints

medianlat=np.median(lats) medianlon=np.median(lons)

maxlat=np.max(lats) maxlon=np.max(lons) minlat=np.min(lats) minlon=np.min(lons) cutlat=(maxlat+minlat)/2 cutlon=(minlat+minlon)/2

global deltarectlon, deltarectlat

rectlon1=minlon1-deltarectlon rectlat1=minlat1-deltarectlat rectlon2=maxlon1+deltarectlon rectlat2=maxlat1+deltarectlat

rectlonc=(rectlon1+rectlon2)/2 rectlatc=(rectlat1+rectlat2)/2

dimlat=rectlat2-rectlat1 dimlon=rectlon2-rectlon1

deltalat=dimy/dimlat deltalon=dimx/dimlon

global suspect global suspectlat, suspectlon, suspectx, suspecty

plon=suspectlon plat=suspectlat

px=int((plon-rectlon1)*deltalon) py=int((plat-rectlat1)*deltalat)

suspectx=px suspecty=py

global data, datax, datay

for n in range(0, numpoints): px=int((lons[n]-rectlon1)*deltalon) py=int((lats[n]-rectlat1)*deltalat) #print(px,py) data.append([px,py]) datax.append(px) datay.append(py)

data=np.array(data) datax=np.array(datax) datay=np.array(datay)

medianx=np.median(datax) mediany=np.median(datay)

#print(data) kenterx=0 kentery=0 for n in range(0, numpoints): kenterx=kenterx+datax[n] kentery=kentery+datay[n]

kenterx=kenterx/numpoints kentery=kentery/numpoints

meanradius=0 rx=0 ry=0

for n in range(0, numpoints): rrx=abs(datax[n]-kenterx) rry=abs(datay[n]-kentery) rx=rx+rrx ry=ry+rry

rx=rx/numpoints ry=ry/numpoints

meanradius=int(math.sqrt(rx*rx+ry*ry))

rx=0 ry=0 rr=0

for n in range(0, numpoints): rx=abs(datax[n]-kenterx) ry=abs(datay[n]-kentery) rra=int(math.sqrt(rx*rx+ry*ry)) rrb=abs(rra-meanradius) rr=rr+rrb

devradius=rr/numpoints

rr=0

global meanpointsdistance global meanpointsdistance_manhattan

mux=0 muy=0 rrm=0

for m in range(0, numpoints): for n in range(0, numpoints): if(n!=m): rx=abs(datax[n]-datax[m]) ry=abs(datay[n]-datay[m]) mux=mux+rx muy=muy+ry rra=int(math.sqrt(rx*rx+ry*ry)) rr=rr+rra rrm=rrm+mux+muy

meanpointsdistance=rr/(((numpoints*numpoints)-numpoints)) meanpointsdistance_manhattan=rrm/(((numpoints*numpoints)-numpoints))

global minimumdistance

minimumdistance=123456789

for m in range(0, numpoints): for n in range(0, numpoints): if(n!=m): rx=abs(datax[n]-datax[m]) ry=abs(datay[n]-datay[m]) rr=int(math.sqrt(rx*rx+ry*ry)) if(rr>0): if(rr<minimumdistance): minimumdistance=rr

#minimumdistance=rr/((numpoints*numpoints))

return(0)

def save_map(map1): global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

delay=5

sw1=[ rectlat1, rectlon1] ne1=[ rectlat2, rectlon2]

map1.fit_bounds([sw1, ne1])

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

map1.save(fn)

#quit(-1) browser = webdriver.Firefox() browser.get(tmpurl) time.sleep(delay) browser.save_screenshot('testmap.png')

browser.quit() return(0)

def create_contours_for_folium(razz1):

global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

geojson1=[]

dimlon=rectlon2-rectlon1 dimlat=rectlat2-rectlat1

deltalon=dimlon/dimx deltalat=dimlat/dimy

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

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

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

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] #ros_mesh = sp.ndimage.filters.gaussian_filter(ros_mesh, sigma, mode='constant') levels=18 #colors = ['blue','royalblue', 'navy','pink', 'mediumpurple', 'darkorchid', 'plum', 'm', 'mediumvioletred', 'palevioletred', 'crimson', # 'magenta','pink','red','yellow','orange', 'brown','green', 'darkgreen'] #colors = ['blue','green','yellow','orange','red', "purple"] colors = ['blue','green','yellow','orange','#3f0000'] levels = len(colors) cm = branca.colormap.LinearColormap(colors, vmin=vmin, vmax=vmax).to_step(levels) #levels=10 #levels=[1,10,100,200] #levels=32 #levels=[0.95,1.0] #lev50=int(255*0.50) #lev75=int(255*0.75) lev90=int(255*0.90) lev95=int(255*0.95) lev97=int(255*0.97) lev98=int(255*0.98) lev100=int(255*1.0) #levels=[lev50,lev75,lev90,lev95,lev97,lev100] levels=[lev90,lev95,lev97,lev98,lev100] #print() contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.5, colors=colors, vmin=vmin, vmax=vmax) geojson1 = geojsoncontour.contourf_to_geojson(contourf=contourf,min_angle_deg=0.01,ndigits=5,stroke_width=3,fill_opacity=0.5)

return(geojson1)

def draw_folium_map(raz1):

global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius global suspectlat, suspectlon global data, lats, lons global minlat, maxlat, minlon, maxlon, cutlon, cutlat global zoomfac1 global drawmatrix global drawimage global imagename1 global imglonsw1,imglatsw1,imglonne1,imglatne1

print("Folium map ...") #plt.imshow(raz1)

#plt.show()

#quit(-1)

sw1=[ rectlat1, rectlon1] ne1=[ rectlat2, rectlon2]

mini1=np.min(raz1) maxi1=np.max(raz1) delti1=maxi1-mini1

raz2=(raz1-mini1)/delti1 raz2=np.exp(raz2) mini1=np.min(raz2) maxi1=np.max(raz2) delti1=maxi1-mini1

raz3=(raz2-mini1)/delti1

raz3=raz3*255

geojson1=create_contours_for_folium(raz3)

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

imr=imr*1 imb=imb*0 img=img*0 ima=imr*0.5 #raz4=np.dstack((imr, img, imb, ima)) raz4=np.dstack((imr, img, imb, ima)) raz5= (raz4* (255.0/raz4.max())).astype(np.uint8)

map1 = folium.Map(location=[kenterlat, kenterlon], TileProvider='OpenStreetMap', zoom_start=zoomfac1, control_scale = True,)

# Plot the contour on Folium map #folium.GeoJson(geojson1,style_function=lambda x: {'color':"#00000",'weight':5.0,'fillColor': "#ff0000",'opacity': 0.5,}).add_to(map1) #folium.GeoJson( geojson1, ).add_to(map1) folium.GeoJson( geojson1, style_function=lambda x: { 'color': x['properties']['stroke'], 'weight': x['properties']['stroke-width']*0.25, 'fillColor': x['properties']['fill'], 'opacity': 0.33, }).add_to(map1)

if(drawmatrix==1): folium.raster_layers.ImageOverlay(raz5,[[ rectlat1, rectlon1],[ rectlat2, rectlon2]],opacity=1,).add_to(map1) folium.LayerControl().add_to(map1)

if(drawimage==1): print("drawimage") # square marker star maybe also icon_square = folium.plugins.BeautifyIcon(icon_shape='rectangle-dot', border_color='red', border_width=10,) #folium.Marker([imglatsw1,imglonsw1 ],).add_to(map1) #folium.Marker([imglatne1, imglonne1] ,).add_to(map1) ## from wiki ... ###imagename1="ripper1.png"

folium.raster_layers.ImageOverlay(imagename1,[[imglatsw1,imglonsw1 ],[imglatne1, imglonne1]],opacity=0.5,).add_to(map1) folium.LayerControl().add_to(map1)

global draw_suspect, draw_skala

lenu=len(lats) for n in range(0, lenu): latx1=lats[n] lonx1=lons[n] folium.Circle([latx1, lonx1], 200*draw_skala, fill_color="blue", opacity=1.0,color = 'blue', fill=True).add_child(folium.Popup('+')).add_to(map1)

latx1=kenterlat lonx1=kenterlon

## NOT OK !!! draw maximum point drawmax=0 ## NOT OK!!! if(drawmax==1): ## !!! NOT OK xlons1=np.linspace(rectlon1,rectlon2,dimx) xlats1=np.linspace(rectlat1,rectlat2,dimy)

maxi1=np.max(raz5) indexes0=np.argwhere(raz5 == maxi1) indexes1=np.array(indexes0[0]) #print("I ", indexes1) #print(indexes1[0]) #print(indexes1[1]) ix1=indexes1[0] iy1=dimy-indexes1[1]

#rectlat1, rectlon1, rectlat2, rectlon2

baselat1=rectlat1 baselon1=rectlon1

dlat1=rectlat2-rectlat1 dlon1=rectlon2-rectlon1

lonco1=dlon1/dimx latco1=dlat1/dimy

#malon1=baselon1+lonco1*ix1 #malat1=baselat1+latco1*iy1 malon1=xlons1[ix1] malat1=xlats1[iy1]

#print("map rect ", rectlat1, rectlon1, rectlat2, rectlon2) print("max point lonlat ", malon1, malat1)

folium.Circle([malat1, malon1], 500*draw_skala, fill_color="green", opacity=1.0,color = 'green', fill=True).add_child(folium.Popup('Williams')).add_to(map1) folium.map.Marker( [malat1 + 0.0001, malon1 - 0.0001], icon=folium.DivIcon( icon_size=(150,36), icon_anchor=(0,0),

) ).add_to(map1)

if(draw_suspect==1):

latx1=suspectlat lonx1=suspectlon folium.Circle([latx1, lonx1], 300*draw_skala, fill_color="black", opacity=1.0,color = 'black', fill=True).add_child(folium.Popup('Williams')).add_to(map1) folium.map.Marker( [latx1 + 0.0001, lonx1 - 0.0001], icon=folium.DivIcon( icon_size=(150,36), icon_anchor=(0,0),

) ).add_to(map1)

#print(" Kenter ", lonx1, latx1)

#folium.Circle([latx1, lonx1], 200*draw_skala, fill_color="green", opacity=1.0,color = 'green', fill=True).add_child(folium.Popup('Barycenter')).add_to(map1) #folium.map.Marker(

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

   #    )
   #).add_to(map1)	

latx1=medianlat lonx1=medianlon

#folium.Circle([latx1, lonx1], 200*draw_skala, fill_color="yellow", opacity=1.0,color = 'yellow', fill=True).add_child(folium.Popup('Median')).add_to(map1) #folium.map.Marker(

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

   #    )
   #).add_to(map1)	

latx1=medianlat lonx1=medianlon

#folium.Circle([latx1, lonx1], 450*draw_skala, fill_color="blue", opacity=1.0,color = 'blue', fill=True).add_child(folium.Popup('MinMax')).add_to(map1) #folium.map.Marker(

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

   #    )
   #).add_to(map1)	
   	

save_map(map1)

return(0)

def view_proba_grid(mata1, L): ## plt ! maybe NOK !!! global rectlat1, rectlat2, rectlon1, rectlon2 global lons, lats global datax, datay

print("View grid ...")

mata1=np.flipud(mata1) #mata1=np.rot90(mata1) #mata1=np.rot90(mata1)

mini1=np.min(mata1) maxi1=np.max(mata1) delti1=maxi1-mini1

mata1=(mata1-mini1)/delti1

fig = plt.figure(dpi=150)

numticks=3 xdim1=np.linspace(0,dimx,numticks) xdim1=np.round(xdim1*100/100) ydim1=np.linspace(0,dimy,numticks) #print(rectlon1) xlabels1=np.linspace(rectlon2, rectlon1, numticks) xlabels1=np.round(xlabels1,4) #print(xlabels1) ylabels1=np.linspace(rectlat1, rectlat2, numticks) ylabels1=np.round(ylabels1,4) ext1=[10,20,80,80] plt.title("Locatiors of cimes - Atlanta serial killer 1979-1981") plt.suptitle("Geographic profiling")

img = plt.imshow(mata1, interpolation='nearest', origin='lower' ) #plt.scatter([a[1] for a in L], [a[0] for a in L], color="w", label ='Crimes', s = 70 ) plt.scatter(datax, datay, color="w", label ='Crimes', s = 45) #plt.contour(proba_matrix, levels=[], color="orange") level1=np.max(np.array(mata1))*0.975 plt.contour(mata1, levels=[level1], colors=["#7f0000"], linewidths=[1], alpha=0.75 ) level1=np.max(np.array(mata1))*0.95 plt.contour(mata1, levels=[level1], colors=["red"], linewidths=[1], alpha=0.75 ) level1=np.max(np.array(mata1))*0.90 plt.contour(mata1, levels=[level1], colors=["orange"], linewidths=[0.5], alpha=0.75 )

global draw_suspect

if(draw_suspect==1): plt.scatter(suspectx, suspecty, color="r", marker="+", label='Williams', s = 500 )

plt.scatter(suspectx, suspecty, color="r", marker="o", s = 60 ) plt.scatter(kenterx, kentery, color="g", marker="x", label='Barycenter', s = 300 )

drawmax=1

if(drawmax==1): maxi1=np.max(mata1) indexes0=np.argwhere(mata1 == maxi1) indexes1=np.array(indexes0[0]) print("I ", indexes1) print(indexes1[0]) print(indexes1[1]) ix1=indexes1[0] iy1=indexes1[1] plt.scatter(iy1,ix1, color="#7f007f", marker="o",facecolors='none', lw=3, alpha=0.95, s=400, label="Maximum probability")

ax = plt.gca() plt.xticks(xdim1, xlabels1, fontsize=9) plt.yticks(ydim1, ylabels1, fontsize=9)

plt.legend(numpoints=1) plt.show() return(0)

1. 1. distance based estimation

def distance_formula(x, y, c, w):

   """
   x: clime location coords
   y: result grid points
   c: decay constant
   w: weights of locs
   """
   dists = cdist(x[:, :2], y, metric='euclidean') #3 ok
   #dists = cdist(x[:, :2], y, metric='chebyshev')  ## ok      
   ##dists = cdist(x[:, :2], y, metric='cityblock')
   #dists = cdist(x[:, :2], y, metric='mahalanobis')

   #dists = cdist(x[:, :2], y, metric='minkowski') #
   #dists = cdist(x[:, :2], y, metric='matching') #   
   weights = w[:, np.newaxis] / (1 + c*dists)
   return np.sum(weights, axis=0)

def calculate_weights(): global dimx, datax, datay, kenterx, kentery weights=[] for n in range(0, numpoints): rx=datax[n]-kenterx ry=datay[n]-kentery rr=math.sqrt(rx*rx+ry*ry) rr=rr/dimx we=1 if(rr>0): we=1-rr weights.append(we)

return(np.array(weights))

def distance_matrix(m,n,L): global dimx, dimy L2=np.array(L) #wei=calculate_weights() #grid_points = np.mgrid[0:dimx:int(dimx), 0:200:200j].reshape(2,-1).T grid_points = np.mgrid[0:dimx, 0:dimy].reshape(2,-1).T ## loc gridpoints distance_values = distance_formula(L2, grid_points, c=0.025, w=np.ones(L2.shape[0])) #distance_values = distance_formula(L2, grid_points, c=0.025, w=wei) # Reshape the Rossmo values into a grid

mata1= distance_values.reshape(dimx, dimy) #mata1=np.flipud(np.rot90(mata1)) mata1=np.fliplr(mata1) mata1=np.rot90(mata1) mata1=np.flipud(mata1) return(mata1)

def normal_matrix(m, n, L): global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

#numpoints=len(lons)

mata1=np.zeros((m, n))

mr=0

for n in range(0, numpoints): ax=abs(lats[n]-kenterx) ay=abs(lons[n]-kentery) ar=int(math.sqrt(ax*ax+ay*ay)) mra=abs(ar-meanradius) mr=mr+mra

## deviation from meanradius mr=mr/numpoints

#deviation=mr ## deviation fro #meanvalue=meanradius

for iy in range(0, dimy): for ix in range(0, dimx): rx=ix-kenterx ry=iy-kentery rr=math.sqrt(rx*rx+ry*ry) mata1[iy,ix]=rr #mata1[iy,ix]=norm.pdf(rr, loc=0, scale=meanradius)

mata1=norm.pdf(mata1, loc=0, scale=meanradius) #mata1=np.fliplr(mata1) #mata1=np.rot90(mata1) mata1=np.flipud(mata1)

return(mata1)

def normal_distribution_ring(sites, buffer_dist, param_f, param_g):

global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

## normal distribution meanr 1*mean radius stdev 5 mata1=np.zeros((dimy, dimx))

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 ## normal, ring radius 1, p0=st.norm.pdf(distance,param_f,param_g) #p0=math.log(math.pow(distance, -1)) p=p+p0

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

mata1=np.flipud(mata1) return(mata1)

def exponential_ring(sites, buffer_dist, param_g, param_f): global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius mata1=np.zeros((dimy, dimx))

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<=1): p=math.exp(-abs(1-distance)/param_g) if(distance>1): p=math.exp(-abs(1-distance)/param_f) p0=p0+p

mata1[y,x]=p0

mata1=np.flipud(mata1) return (mata1)

def exponential_distance(sites, buffer_dist):

global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

## normal distribution meanr 1*mean radius stdev 5 mata1=np.zeros((dimy, dimx))

p0=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 base1=(math.sqrt(ux*ux+uy*uy))/buffer_dist p=math.exp(-1*base1) p0=p0+p

mata1[y,x]=p0 ## tsek it!

mata1=np.flipud(mata1) return(mata1)

def normalize_matrix(inmat1): min1=np.min(inmat1) max1=np.max(inmat1) del1=max1-min1 outmat1=(inmat1-min1)/del1 return(outmat1)

def density_matrix(): global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius

weights=np.zeros(numpoints)+1

funk1='multiquadric' #funk1='inverse' #funk1='gaussian' #funk1='linear' #funk1='cubic' #funk1='quintic' #funk1='thin_plate'

## eplslon smooth norm mode rbf_adj = Rbf(datax, datay, weights, function=funk1)

x_fine = np.linspace(0, dimx, dimx) y_fine = np.linspace(0, dimy, dimy)

x_grid, y_grid = np.meshgrid(x_fine, y_fine)

mata0 = rbf_adj(x_grid.ravel(), y_grid.ravel()).reshape(x_grid.shape)

mata1=np.array(mata0)

mata1=np.flipud(mata1)

return(mata1)

def save_netcdf(outfilename1, mata1): global numpoints, dimx,dimy,rectlon1,rectlon2,rectlat1,rectlat2,rectlatc,rectlatc global lats, lons, latsize,lonsize,dlat,dlon,kenterlat,kenterlon,kenterx,kentery,meanradius print("Saving netcdf4 ...") ncout2 = netCDF4.Dataset(outfilename1,mode='w',format='NETCDF4_CLASSIC') lat_dim = ncout2.createDimension('lat', dimy) lon_dim = ncout2.createDimension('lon', dimx) ncout2.title='p' ncout2.subtitle="value" lat = ncout2.createVariable('lat', np.float32, ('lat',)) lat.units = 'degrees_north' lat.long_name = 'latitude' lon = ncout2.createVariable('lon', np.float32, ('lon',)) lon.units = 'degrees_east' lon.long_name = 'longitude' temp = ncout2.createVariable('P',np.float64,('lat','lon')) temp.units = 'value' temp.standard_name = 'P' temp.standard_name = 'Probability' nlats = len(lat_dim) nlons = len(lon_dim) temp[:,:] = np.flipud(mata1) latdim1=rectlat2-rectlat1 londim1=rectlon2-rectlon1 latsi1 = rectlat1 + (latdim1/dimy)*np.arange(dimy) lonsi1 = rectlon1 + (londim1/dimx)*np.arange(dimx) #print(lonsi1) lat[:] = latsi1 lon[:] = lonsi1

ncout2.close() return(0)

def calculate_probability(L,A):

   global meanradius, devradius, meanpointsdistance, minimumdistance
   global meanpointsdistance_manhatta
   global bufferi, f, g

   m=dimy
   n=dimx
   
   ## development only
   ## drop crime locations all but certain number near barycenter
   #crime_locations_dropper(6)
   
   #quit(-1)
   
   generate_coords()
   
   print(crimelocations)
   
   bufferi=meanpointsdistance/2
   
   print("")
   print("Mean distance ", meanpointsdistance)
   print("Bufferi       ", bufferi)        
   print("f             ", f) 
   print("g             ", g)
   print("")      
   #quit(-1)
   
   L=list(data)
   
   mata1 = rossmo_matrix(m,n,L) ## ok basic rossmo, seems unaccurate
   #koeff1 = 35
   #mata1=rossmo_2(koeff1,f,g,bufferi, L) ##      
   #mata1=rossmo_3(L, bufferi, f,g)  ##
   #mata1=exponential_ring(L, bufferi, f, g) ## nok      
   #mata1=exponential_ring(L, bufferi, 2, 2) ## nok      
   #mata1=density_matrix()
   
   #mata1=exponential_distance(L, devradius*16) 
      
   #mata1=distance_matrix(m,n,L) ## OK distance from different points
   #mata1=normal_matrix(m, n, L) ## is ok nr not )= distance from barycenter
 
   #mata1=normal_distribution_ring(L,devradius/2, 1,0.25)

   #mata1=normalize_matrix(mata1)  
   
   save_netcdf("test.nc", mata1)
   
   nm1=normalize_matrix(mata1)
   nm2=np.uint8(nm1*255)
   im1 = Image.fromarray(nm2,'L')
   #im1.show()

   im1.save("analysis1.png") 
   
   # maybe NOK ???
   #view_proba_grid(mata1, L)
   
   draw_folium_map(mata1)

   return(0)
   

1. 1. create estimation of propability of location of suspect

calculate_probability(crimelocations, suspect)