File:Rossmo analysis example desalvo murders 1.png

Size of this preview: 800 × 578 pixels. Other resolutions: 320 × 231 pixels | 640 × 463 pixels | 1,047 × 757 pixels.
Original file ‎(1,047 × 757 pixels, file size: 169 KB, MIME type: image/png)
Captions

English
Rossmo analysis example: DeSalvo murders
Summary

Description	English: Gaographic analysis of serial crime: attempt to guess location of residence of killer.
Date	14 September 2023
Source	Own work
Author	Merikanto
Sources of equations and data
1. independly coded, but based formulas on
2. rossmo equations from
3. https://medium.com/@errazkim/how-to-catch-serial-killers-using-one-simple-formula-29b1cbd58963
https://sites.math.washington.edu/~morrow/mcm/7272.pdf
Why Crime Doesn’t Pay: Locating Criminals Through Geographic Profiling
Control Number: #7272
February 22, 2010
1. richard chase data from
1. https://jeremykun.com/2011/07/20/serial-killers/
2. Hunting Serial Killers “Tonight’s the Night”
Posted on July 20, 2011 by Jeremy Kun
Python3 source code
#########################
#
## attempt to make rossmo and gaussian rossmooth crime location analysis
## python3 code
##
# 20.9.2023 0000.0011
#
#############################

import os, time
import numpy as np
import math
from math import sqrt, exp

import scipy.stats as stats
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
import matplotlib.pyplot as plt
from matplotlib.ticker import StrMethodFormatter

from matplotlib import ticker

## image and raster save routines
from PIL import Image
import netCDF4

## map routines

import folium
from folium import plugins
from selenium import webdriver
import branca
import geojsoncontour



#
## independly coded, but based formulas on
## rossmo equations from
## https://medium.com/@errazkim/how-to-catch-serial-killers-using-one-simple-formula-29b1cbd58963
#
# https://sites.math.washington.edu/~morrow/mcm/7272.pdf
# Why Crime Doesn’t Pay: Locating Criminals Through Geographic Profiling
#Control Number: #7272
#February 22, 2010
#
## richard chase data from
#
## https://jeremykun.com/2011/07/20/serial-killers/
## Hunting Serial Killers “Tonight’s the Night”
# Posted on July 20, 2011 by Jeremy Kun	
#
##


buf1=0 ## automatically adjusted

#f=.334 ## often 0.5
#f=.333 ## often 0.5
#f=0.5
#f=1.2
#g=.667 ## 0.1 ... 1.25

#g=0.67 ## 0.1 ... 1.25

g=1
f=1
## f 0.5 g 0.67

#f=0.5 ## nok

#g=0.66

#f=1
#g=1


koo=1


sizex=256
sizey=256

#sizex=100
#sizey=100


#coordsx=np.array([50,20,80,60])*1
#coordsy=np.array([20,60,10,80])*1

#sites1=26
#coordsx=np.random.rand(sites1)*sizex
#coordsy=np.random.rand(sites1)*sizey

#coordsx=np.array([0+20,100-20,100-20,0+20])*2
#coordsy=np.array([0+20,0+20,100-20,100-20])*2

#coordsx=np.array([50,50,1,99])*2
#coordsy=np.array([1,99,50,50])*2


#residencex1=50*2+1
#residencey1=50*2+1
#residencex2=50*2-1
#residencey2=50*2-1


   
richardchase_coordsx = np.array([3, 15, 19, 21, 25 ])*2  
richardchase_coordsy = np.array([17, 3, 27, 22, 18 ])*2
   
richardchase_residencex=19*2
richardchase_residencey=17*2


desalvo_coordsx = np.array([10,13, 15, 17,18,18, 19, 19, 20, 20, 20, 29, 33 ])

desalvo_coordsy =np.array([ 48, 8, 11 , 8 , 7  , 9,  4,  8, 9, 10, 11,23,28])

desalvo_residencex = 19
desalvo_residencey = 18



sutcliffe_coordsx =np.array([ 5,8, 50,53,56, 59,62, 63, 63,  64,  69, 73,  80,  81,  83,  83,  85, 85, 90])

sutcliffe_coordsy = np.array([ 1, 7, 99,  68,  72, 59,  57,  85,  87,  83,  82,88,  88,  89,  88,  87,  85, 83,90 ])


sutcliffe_residencesx = np.array([60, 58])
sutcliffe_residencesy = np.array([88, 81])



coordsx=desalvo_coordsx*4 
coordsy=desalvo_coordsy*4

residencex1=desalvo_residencex*4
residencey1=desalvo_residencey*4
residencex2=desalvo_residencex*4
residencey2=desalvo_residencey*4
   
   
   
#coordsx=sutcliffe_coordsx*2 
#coordsy=sutcliffe_coordsy*2


#residencex1=sutcliffe_residencesx[0]*2
#residencey1=sutcliffe_residencesy[0]*2
#residencex2=sutcliffe_residencesx[1]*2
#residencey2=sutcliffe_residencesy[1]*2


#coordsx=richardchase_coordsx*4 
#coordsy=richardchase_coordsy*4 


#residencex1=richardchase_residencex*4
#residencey1=richardchase_residencey*4
#residencex2=richardchase_residencex*4
#residencey2=richardchase_residencey*4   

#
## locations from
# https://medium.com/@errazkim/how-to-catch-serial-killers-using-one-simple-formula-29b1cbd58963
# Catching serial killers using Rossmo’s formula 
# Understanding Rossmo’s formula with Atlanta’s murders
# Mohamed Errazki
# Jul 1, 2022
#

#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]])

## suspect locations, available in some cases
## Atlanta Williams

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

#extent1<-c(-84.946577,-83.946577,33.254065,34.254065)
 

## atlanta murders 1079-1981 williams
atlanta_coordsx=np.array([106,115,144,127,92,152,115,123,152,137,121,132,150,125,95,225,157,178,209,78,67,134,68,114,158,114
])
atlanta_coordsy=np.array([114,103,124,116,114,126,104,128,112,124,141,108,108,154,109,109,149,109,96,110,102,131,102,140,124,140
])

atlanta_residencex=127
atlanta_residencey=127

#coordsx=atlanta_coordsx*1 
#coordsy=256-atlanta_coordsy*1


#residencex1=atlanta_residencex*1
#residencey1=atlanta_residencey*1
#residencex2=atlanta_residencex*1
#residencey2=atlanta_residencey*1 

#numpoints=len(coordsx)
numpoints=0


centerx=0
centery=0
deviationx=0
deviationy=0

mindist=0

standard_distance_manhattan=0
standard_distance_euclidian=0



## jact the ripper, beta only

#jack_the_ripper_lats=np.array([51.5203785, 51.5186270, 51.5199207, 51.5137500,51.5137633])
#jack_the_ripper_lons=np.array([-0.0737000,-0.0752304,-0.0610239,-0.0779247,-0.0657004])


#lons=jack_the_ripper_lons
#lats=jack_the_ripper_lats

atlanta_lats=np.array([33.703093,33.660032,33.741141,33.711061,33.701493,33.746652,33.6605585,33.755227,33.692281,33.738875,33.805397,33.677783,33.679033,33.858629,33.68205,33.68164,33.837342,33.680747,33.631259,33.683782,33.653495,33.766465,33.653852,33.802901,33.7392226,33.804113])
 

atlanta_lons = np.array([-84.532406,-84.49509,-84.383959,-84.447227,-84.584169,-84.350482,-84.4941276,-84.465294,-84.350066,-84.408613, -84.470401,-84.427292,-84.358048, -84.455166,-84.573247,-84.067554,-84.332364,-84.249387,-84.128966,-84.64159,-84.681008,-84.422132,-84.6803,-84.500141,-84.3287373,-84.499154 ])



atlanta_suspect_lat=33.754065
atlanta_suspect_lon=-84.446577

numpoints0=0

lons0=atlanta_lons
lats0=atlanta_lats

lons=atlanta_lons
lats=atlanta_lats


suspectlat=atlanta_suspect_lat
suspectlon=atlanta_suspect_lon



centerlat=1
centerlon=1
crimeboxlon1=1
crimeboxlat1=1
crimeboxlon2=1
crimeboxlat2=1


#coordsx=[] 
#coordsy=[]


## output analysis raster or matrix

raster=np.zeros((sizey, sizex))

def crime_locations_dropper(selekt1):

	global numpoints,numpoints0, lons, lats


	lons0=np.copy(lons)
	lats0=np.copy(lats)	
	numpoints0=len(lons)-1
	#print(numpoints0)	
	
	cenx=0
	ceny=0
	rx=0
	ry=0
	
	for n in range(0, numpoints0):
		rax=lons0[0]
		ray=lats0[0]
		rx=rx+rax
		ry=ry+ray
	
	cenx=rx/numpoints0
	ceny=ry/numpoints0
	
	#print(cenx, ceny)
	
	
	rtab0=[]
	utab0=[]
	
	for n in range(0, numpoints0):
		x=lons0[n]
		y=lats0[n]
		rax=x-cenx
		ray=y-ceny
		rr=math.sqrt(rax*rax+ray*ray)
		#print(rr, rax, 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)
	
	#quit(-1)
				
	
	#print(utab3)	
	utab4=utab3[:,1:].T
	#print(utab4)
	
	lons=None
	lats=None
	
	lons=utab4[0]
	lats=utab4[1]
	
	
	#print(lons)
	#print(lats)
	#crimelocations=None
	
	#crimelocations=utab4		
	#koordis=np.copy(crimelocations)
	
	return(0)



def density_matrix_1():
	global numpoints
	global sizex, sizey,coordsx, coordsy, centerx, centery
	global sizex, sizey, raster
	global crimeboxlon1,crimeboxlat1,crimeboxlon2,crimeboxlat2
	
	weights=np.zeros(numpoints)+1
	
	ifunk1='multiquadric'
	#funk1='inverse'
	#funk1='gaussian'
	#funk1='linear'
	#funk1='cubic'
	#funk1='quintic'
	#funk1='thin_plate'

	## eplslon smooth norm mode
	rbf_adj = Rbf(coordsx, coordsy, weights, function=ifunk1)

	x_fine = np.linspace(0, sizex, sizex)
	y_fine = np.linspace(0, sizey, sizey)

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

	raster = rbf_adj(x_grid.ravel(), y_grid.ravel()).reshape(x_grid.shape)
	
	#mata1=np.array(mata0)
	
	#mata1=np.flipud(mata1)
	
	return(0)








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

    #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 sizex, sizey,coordsx, coordsy, centerx, centery
 
	weights=[]
	for n in range(0, numpoints):
		rx=coordsx[n]-centerx	
		ry=coordsy[n]-centery
		rr=math.sqrt(rx*rx+ry*ry)
		rr=rr/sizex
		we=1
		if(rr>0): we=1-rr
		weights.append(we)
	
	return(np.array(weights))
			


def distance_matrix():
	global coordsx, coordsy,sizex, sizey, raster
	
	wei=calculate_weights()
	#grid_points = np.mgrid[0:dimx:int(dimx), 0:200:200j].reshape(2,-1).T
	grid_points = np.mgrid[0:sizex, 0:sizey].reshape(2,-1).T
	## loc gridpoints
	inpoints=np.array(np.dstack((coordsx, coordsy)))[0]
	#print(np.shape(inpoints))
	##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)
	## set all weights to 1:
	wei=wei*1
	distance_values = distance_formula(inpoints, grid_points, c=0.025, w=wei)	
	# Reshape the Rossmo values into a grid
	
	raster=distance_values.reshape(sizex,sizey)

	return(0)


def save_netcdf(outfilename1):
	## save netcdf: WARNING have not crs!!!
	global sizex, sizey, raster
	global crimeboxlon1,crimeboxlat1,crimeboxlon2,crimeboxlat2

	print("Saving netcdf4 ...")
	ncout2 = netCDF4.Dataset(outfilename1,mode='w',format='NETCDF4_CLASSIC')
	lat_dim = ncout2.createDimension('lat', sizex) 
	lon_dim = ncout2.createDimension('lon', sizey) 
	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(raster)
	latdim1=crimeboxlat2-crimeboxlat1
	londim1=crimeboxlon2-crimeboxlon1	
	latsi1 = crimeboxlat1 + (latdim1/sizex)*np.arange(sizey)  
	lonsi1 = crimeboxlon1 + (londim1/sizey)*np.arange(sizex) 
	#print(lonsi1)
	lat[:] = latsi1  
	lon[:] = lonsi1

	ncout2.close()
	return(0)






## boundary grs method
def calculate_grs_d3(dx,dy,fii, baffer1, cee1):
	dx=dx
	dy=dy
	
	a0=1-fii
	
	a1=abs(dx)+abs(dy)
	
	#print(a1)
	
	a2=(baffer1)-a1
	
	#a2=a2/100
	
	
	#print(a2)
	
	
	a3=-1*a2*a2
	
	
	#print(a3)
	
	a4=exp(a3)
	
	#print("exp (a4)", a4)
	
	a5=cee1*a4
	
	d2=a0*a5
	
	#print("dee2",d2)
	
	return(d2)









## not ok not ok not ok
def gaussian_boundary_method():
	## attempt to make GRS, Gaussian RosSmooth ...
	## warning  NOT OK!!!!
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global centerx,centery,deviationx,deviationy	
	global standard_distance_manhattan
	global standard_distance_euclidian	
	
	
	
	#buf1=standard_distance_manhattan
	#buf1=mindist/2
	buf1=standard_distance_euclidian/2
	
	### "baffer", bot buffer for tesing GRS boundary
	baffer1=mindist
	
	## epsilon 0.5

	epsilon=0.5
	
	
	print("Gaussian boundary ")
	
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	
	## tweaked
	#buf1=mindist
	
		
	
	#quit(-1)
		
	
	
	## grs rossmo  test 
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
					
			
			for n in range(0,numpoints):
				
				px=coordsx[n]
				py=coordsy[n]
				dx=px-ix
				dy=py-iy
				dis=abs(dx)+abs(dy)
				
				##ensure that distance causes not overflow!
				if(dis==0): dis=1
								
				## smoothed rossmo
				#fii=0
				#if(dis>buf1): fii=1
				##fii=epsilon
				#fii=0.5
				
				cee1=1
				fii=0.5
				cap1=fii/pow(epsilon, f)
				
				
				cee1=cap1
				

				#term1=calculate_grs_d1(dx,dy,f, fii, epsilon)
				#fii=0.5
				#term2=calculate_grs_d2(dx,dy,fii, buf1, cee1)
				#gaussian1=gaussian_term(ix,iy,fii,epsilon)						
				term3=calculate_grs_d2(dx,dy,fii, baffer1, cee1)				
				#print(term2)
				#quit(-1)
				
				#koo2=1

				val1=term3
				
				#val1=gaussian1
				#val1=term1+gaussian1
				#val1=term2
								
				p=val1
				
				pp=pp+p
							
			
			raster[iy,ix]=pp
			
					
		
	return(0)













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

def create_contours_for_folium(rasta0):
	
	global numpoints, sizex, sizey
	global centerlon, centerlat
	global lons,lats,suspectlat,suspectlon
	global centerlat,centerlon
	global crimeboxlon1,crimeboxlat1,crimeboxlon2,crimeboxlat2
	
	

	rasta=rasta0
	
	#plt.imshow(rasta0)
	#plt.show()
	#quit(-1)
			
	geojson1=[]
	
	dimlon=crimeboxlon2-crimeboxlon1
	dimlat=crimeboxlat2-crimeboxlat1
	
	
	deltalon=dimlon/sizex
	deltalat=dimlat/sizey
	
	lontab=[]
	lattab=[]
	rostab=[]

	for iy in range(0, sizey-1):
		tlat=crimeboxlat2-deltalat*iy	
		for ix in range(0, sizex-1):
			tlon=crimeboxlon1+deltalon*ix
			tros=rasta[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), sizex)
	lat_arr          = np.linspace(np.min(lattab), np.max(lattab), sizey)
	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','#100000']
	
	#levels=20
	
	levels = len(colors)
	#cm     = branca.colormap.LinearColormap(colors, vmin=vmin, vmax=vmax).to_step(levels)
	#levels=10
	lev50=int(255*0.50)
	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]
	#levels=10
	#print()
	
	#plt.imshow(rasta)
	
	
	contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.5, colors=colors, vmin=vmin, vmax=vmax)
	#contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.5, vmin=vmin, vmax=vmax)
	

	
	#plt.show()
	
	
	geojson1 = geojsoncontour.contourf_to_geojson(contourf=contourf,min_angle_deg=0.01,ndigits=5,stroke_width=3,fill_opacity=0.5)


	return(geojson1)



## NOK
def draw_folium_map(rasta0):
	
	global numpoints
	global centerlon, centerlat
	global lons,lats,suspectlat,suspectlon
	global centerlat,centerlon
	global crimeboxlon1,crimeboxlat1,crimeboxlon2,crimeboxlat2	
	
	
	zoomfac=12
	
	


	rasta=np.array(normalize_matrix(rasta0)*255).astype(int)
	
	#plt.imshow(rasta)
	
	#plt.show()
	

	geojson1=create_contours_for_folium(rasta)
	
	#quit(-1)
	
	rast2=np.exp(rasta0)
	
	rast2a=np.copy(rast2)*255
	rast2r=np.copy(rast2)*255	
	rast2g=np.copy(rast2)*255*0		
	rast2b=np.copy(rast2)*255*0
	
	rastafari=np.dstack((rast2r,rast2g, rast2b, rast2a))
			
	
	
	
	
	print("Folium map ...")
	
	tilesfrom1='OpenStreetMap'
	##tilesfrom1='StamenToner' NOK
	#tilesfrom1='https://tiles.stadiamaps.com/tiles/stamen_toner_background/{z}/{x}/{y}{r}.png'
	
	map1 = folium.Map(location=[centerlat, centerlon], TileProvider=tilesfrom1, zoom_start=zoomfac, control_scale = True,)

	draw_contour=1
	if(draw_contour==1):
		
		folium.GeoJson(
			geojson1,
			style_function=lambda x: {
				'color':     x['properties']['stroke'],
				'weight':    x['properties']['stroke-width']*1,
				'fillColor': x['properties']['fill'],
				'opacity':   1,
				}).add_to(map1)
			
	drawmatrix=1
	if(drawmatrix==1):
		rasta2=np.flipud(rasta)*10
		
		folium.raster_layers.ImageOverlay(rastafari,[[crimeboxlat1, crimeboxlon1 ],[ crimeboxlat2, crimeboxlon2 ]],opacity=0.5,).add_to(map1)
		folium.LayerControl().add_to(map1)

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

	draw_suspect=1
	if(draw_suspect==1):
		
		latx1=suspectlat
		lonx1=suspectlon
		folium.Circle([latx1, lonx1], 300*1, 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),
			html='<div style="font-size: 14pt">%s</div>' % "Suspect",
			)
		).add_to(map1)


	draw_center=1
	if(draw_center==1):
		
		latx1=centerlat
		lonx1=centerlon
		folium.Circle([latx1, lonx1], 300*1, fill_color="green", opacity=1.0,color = 'green', 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),
			html='<div style="font-size: 14pt">%s</div>' % "Barycenter",
			)
		).add_to(map1)	
	fn='testmap.html'
	tmpurl='file://{path}/{mapfile}'.format(path=os.getcwd(),mapfile=fn)
	
	map1.save(fn)
	
	#quit(-1)
	browser = webdriver.Firefox()
	browser.get(tmpurl)
	delay=5
	time.sleep(delay)
	browser.save_screenshot('testmap.png')

	browser.quit()
	
	return(0)

def translate_geographical_coordinates():
	global numpoints, lons, lats, sizex, sizey, suspectlat, suspectlon
	global crimeboxlon1,crimeboxlat1,crimeboxlon2,crimeboxlat2
	global coordsx,coordsy,residencex1,residencey1,residencex2,residencey2
	global centerlon, centerlat

	numpoints=len(lons)
	
	maxlat=np.max(lats)
	minlat=np.min(lats)
	maxlon=np.max(lons)
	minlon=np.min(lons)
	
	dimlat=maxlat-minlat
	dimlon=maxlon-minlon
	
	crimeboxlon1=minlon-(dimlon/4)
	crimeboxlat1=minlat-(dimlat/4)
	crimeboxlon2=maxlon+(dimlon/4)
	crimeboxlat2=maxlat+(dimlat/4)
	
	centerlon=(crimeboxlon2+crimeboxlon1)/2
	centerlat=(crimeboxlat2+crimeboxlat1)/2
		
	dimboxlon=crimeboxlon2-crimeboxlon1
	dimboxlat=crimeboxlat2-crimeboxlat1
	
	loncoeff=sizex/dimboxlon
	latcoeff=sizey/dimboxlat
		
	for n in range(0, numpoints):
		dlat1=lats[n]-crimeboxlat1
		dlon1=lons[n]-crimeboxlon1
		rix=int(dlon1*loncoeff)
		riy=int(dlat1*latcoeff)	
		print(n,rix,riy)
		coordsx.append(rix)
		coordsy.append(riy)
	
	
	coordsx=np.array(coordsx)
	coordsy=np.array(coordsy)	
	
	residencex1=sizex/2
	residencey1=sizey/2

	
	residencex1=int(loncoeff*(suspectlon-crimeboxlon1))
	residencey1=int(latcoeff*(suspectlat-crimeboxlat1))
	
	print(" Residence ", residencex1, residencey1)

	residencex2=residencex1
	residencey2=residencey1
	
	return(0)



def gaussian_distance_from_center_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global standard_distance_manhattan
	global standard_distance_euclidian
	global centerx, centery	
		
	disk=mindist/2
	print(disk)
	## official
	#buf1=disk
	
	#buf1=standard_distance_manhattan
	buf1=standard_distance_euclidian
	#buf1=buf1/sqrt(2)
	#buf1=100
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	## tweaked
	#buf1=mindist
	
	
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
			#for n in range(0,numpoints):
			px=centerx
			py=centery
			dx=abs(px-ix)
			dy=abs(py-iy)
			#dis=dx+dy
			dis=sqrt(dx*dx+dy*dy)
			disdelta=buf1-dis
			absdisdelta=abs(disdelta)
			
			#discoeff=dis/buf1
			discoeff1=dis/buf1
			
			absdiscoeff=abs(discoeff1)
			gee=1/2/buf1
			fff=2/3/buf1
			val1=0
			#val1=exp(-absdisdelta*gee)
			
			## normal dist
			a1=stats.norm.pdf(discoeff1, loc=0, scale=1)
			val1=a1

			pp=val1
				
			#pp=pp+p
										
			raster[iy,ix]=pp
			
					
		
	return(0)




def exponential_distance_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global standard_distance_manhattan
	global standard_distance_euclidian
	global centerx, centery	
		
	disk=mindist/2
	print(disk)
	## official
	#buf1=disk
	
	#buf1=standard_distance_manhattan
	buf1=standard_distance_euclidian
	#buf1=buf1/sqrt(2)
	#buf1=100
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	## tweaked
	#buf1=mindist
	
	
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
			for n in range(0,numpoints):
				px=coordsx[n]
				py=coordsy[n]
				dx=abs(px-ix)
				dy=abs(py-iy)
				#dis=dx+dy
				dis=sqrt(dx*dx+dy*dy)
				disdelta=buf1-dis
				absdisdelta=abs(disdelta)
				
				discoeff=dis/buf1
				
				absdiscoeff=abs(discoeff)
				gee=1/2/buf1
				fff=2/3/buf1

				val1=0

				a1=exp(-discoeff*1)
				
				## normal dist
				#a1=stats.norm.pdf(discoeff, loc=1, scale=1/2)
				val1=a1
			
				
				p=val1
				
				pp=pp+p
							
			
			raster[iy,ix]=pp
			
					
		
	return(0)










def calculate_grs_d2(dx,dy,fii, buffer1, cee1):
	
	dx=dx
	dy=dy
	
	a0=1-fii
	
	a1=abs(dx)+abs(dy)
	
	#print(a1)
	
	a2=(2*buffer1)-a1
	
	##a2=a2/100
	
	
	#print(a2)
	
	
	a3=-1*a2*a2
	
	
	#print(a3)
	
	a4=exp(a3)
	
	#print("exp (a4)", a4)
	
	a5=cee1*a4
	
	d2=a0*a5
	
	#print("dee2",d2)
	
	return(d2)


def calculate_grs_d1(dx,dy,f, fii, epsilon):
	dx=dx
	dy=dy
	dii=abs(dx)+abs(dy)
	if(dii==0): dii=1
	
	d11=fii/pow(dii,f)
	d12=fii/pow(epsilon, f)
	d1=min(d11,d12)
	
	return(d1)


def gaussian_term(ix,iy,fii,epsilon):
	global centerx,centery,deviationx,deviationy	
	global f
	global xsize
	## gaussian
					
	dex=ix-centerx
	dey=iy-centery
	dex2=dex*dex
	dey2=dey*dey
	bx1=dex2/(2*deviationx)
	by1=dey2/(2*deviationy)
	gex1=-1*(bx1+by1)
	gauss1=exp(gex1) 
							
	## epsilon 0.5
	aaa1=(2*fii)/pow(epsilon,f)
				
	gauss2=aaa1*gauss1
	
	return(gauss2)


def gaussian_rossmooth_1():
	## attempt to make GRS, Gaussian RosSmooth ...
	## warning maybe NOT OK!!!!
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global centerx,centery,deviationx,deviationy	
	global standard_distance_manhattan
	global standard_distance_euclidian	
	
	
	
	#buf1=standard_distance_manhattan
	buf1=mindist
	#buf1=standard_distance_euclidian/2
	
	## epsilon 0.5

	epsilon=0.5
	
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	
	## tweaked
	#buf1=mindist
	
		
	
	#quit(-1)
		
	
	
	## grs rossmo  test 
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0

			## fii=0.5		
			px1=coordsx[0]
			py1=coordsy[0]
			dx1=abs(px1-ix)
			dy1=abs(py1-iy)
			dis=dx1+dy1
			if(dis==0): dis=1			
			
			## smoothed rossmo


			
			#cee1=1
			w1=1
			fii=0
			if(dis>buf1): fii=1
			
			##fii=epsilon
			fii=0.5
			
			
			cap1=fii/pow(epsilon, f)
			cee1=cap1				
			
			dw1=calculate_grs_d1(dx1,dy1,f, fii, epsilon)
			dw2=calculate_grs_d2(dx1,dy1,fii, buf1, cee1)
				
			dewe1=w1*(dw1+dw2)
			
			
			gaussian1=gaussian_term(ix,iy,fii,epsilon)			
			
			
			for n in range(1,numpoints):
				
				px=coordsx[n]
				py=coordsy[n]
				dx=abs(px-ix)
				dy=abs(py-iy)
				dis=dx+dy
				
				##ensure that distance causes not overflow!
				if(dis==0): dis=1
								
				## smoothed rossmo
				fii=0
				if(dis>buf1): fii=1
				##fii=epsilon
				#fii=0.5
				
				cee1=1
				
				cap1=fii/pow(epsilon, f)
				
				
				cee1=cap1
				
				fii=0.5
				term1=calculate_grs_d1(dx,dy,f, fii, epsilon)
				fii=0.5
				term2=calculate_grs_d2(dx,dy,fii, buf1, cee1)					
				
				#print(term2)
				#quit(-1)
				
				#koo2=1

				val1=term1+term2+dewe1+gaussian1
				
				#val1=gaussian1
				#val1=term1+gaussian1
				#val1=term2
								
				p=val1
				
				pp=pp+p
							
			
			raster[iy,ix]=pp
			
					
		
	return(0)




def exponential_ring_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global standard_distance_manhattan
	global standard_distance_euclidian	
		
	disk=mindist/2
	print(disk)
	## official
	#buf1=disk
	
	#buf1=standard_distance_manhattan
	buf1=standard_distance_euclidian
	buf1=buf1/sqrt(2)
	#buf1=100
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	## tweaked
	#buf1=mindist
	
	
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
			for n in range(0,numpoints):
				px=coordsx[n]
				py=coordsy[n]
				dx=abs(px-ix)
				dy=abs(py-iy)
				#dis=dx+dy
				dis=sqrt(dx*dx+dy*dy)
				disdelta=buf1-dis
				absdisdelta=abs(disdelta)
				
				buf1=3
				
				discoeff=dis/buf1
				
				absdiscoeff=abs(discoeff)
				gee=1/2/buf1
				fff=2/3/buf1

				val1=0

				val1=exp(-absdisdelta*gee)
				
	
				
				p=val1
				
				pp=pp+p
							
			
			raster[iy,ix]=pp
			
					
		
	return(0)


def gaussian_ring_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global standard_distance_manhattan
	global standard_distance_euclidian	
		
	disk=mindist/2
	print(disk)
	## official
	#buf1=disk
	
	#buf1=standard_distance_manhattan
	buf1=standard_distance_euclidian
	buf1=buf1/sqrt(2)
	#buf1=100
	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	## tweaked
	#buf1=mindist
	
	
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
			for n in range(0,numpoints):
				px=coordsx[n]
				py=coordsy[n]
				dx=abs(px-ix)
				dy=abs(py-iy)
				#dis=dx+dy
				dis=sqrt(dx*dx+dy*dy)
				disdelta=buf1-dis
				absdisdelta=abs(disdelta)
				
				discoeff=dis/buf1
				if(discoeff==0): discoeff=0.001
				
				absdiscoeff=abs(discoeff)
				gee=1/2/buf1
				fff=2/3/buf1

				val1=0

				#val1=exp(-absdisdelta*gee)
				
				## normal dist
				#a1=stats.norm.pdf(discoeff, loc=1, scale=1/2)
				mu=1
				sigma=1/3
				
				h1=math.pow(abs(discoeff-mu),2)/(2*sigma*sigma)
				h2=math.exp(-1*h1)
				h3=sigma*math.sqrt(2*math.pi)
				
				a1=h2/h3
				
				
				#a1=math.exp(-(discoeff*discoeff)/2)/math.sqrt(2*math.pi)
				
				
				
				
				
				val1=a1
			
				
				p=val1
				
				pp=pp+p
							
			
			raster[iy,ix]=pp
			
					
		
	return(0)


def find_nearest():
	global coordsx,coordsy
	global centerx,centery,deviationx,deviationy	



	
	dis1=9999999
	
	
	for m in range(0,numpoints):
		px=coordsx[m]
		py=coordsy[m]
		for n in range(0,numpoints):
			qx=coordsx[n]
			qy=coordsy[n]
			dx=qx-px
			dy=qy-py
			if(n!=m):
				## manhattan distance
				d1=abs(dx)+abs(dy)
				if(d1>0):
					if(d1<dis1): dis1=d1
				
	px=0
	py=0
	
	for n in range(0,numpoints):
		px=px+coordsx[n]
		py=py+coordsy[n]
	
	centerx=px/numpoints
	centery=py/numpoints
	
	px=0
	py=0
	
	for n in range(0,numpoints):
		px=px+abs(coordsx[n]-centery)
		py=py+abs(coordsy[n]-centerx)	
	
	deviationx=px/numpoints
	deviationy=py/numpoints
	
	
	global standard_distance_manhattan
	global standard_distance_euclidian
		
	standard_distance_manhattan=abs(deviationx)+abs(deviationy)		
	standard_distance_euclidian=sqrt(deviationx*deviationx+deviationy*deviationy)	
	
	#print("center       ", centerx, centery)
	#print("deviation    ",deviationx, deviationy)
	#print("min distance ", dis1)
	
	return(dis1)


def draw_raster_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global residencex, residencey
	global f, g, buffer1
	
	
	raster1=np.copy(raster)
	
	max1=np.max(raster1)
	min1=np.min(raster1)
	del1=max1-min1
	raster1=(raster1-min1)/del1
	
	colors1=list(["blue", "green", "yellow","orange", "red", "violet"])
	#colors1=colors1.reverse()
	plt.suptitle("Rossmo analysis", size=18)
	sup1="Buffer="+str(int(buf1))+ "  f=" +str(round(f,3)) +"  g="+ str(round(g,3))
	plt.title(sup1, size=16)
	
	#plt.rcParams.update({'font.size': 16})

	
	#cmap = plt.get_cmap('gist_ncar')
	
	#plt.set_cmap(cmap)
	#plt.imshow(raster1,origin="lower",  cmap=matplotlib.colormaps['summer'])
	plt.imshow(raster1,origin="lower",  cmap=matplotlib.colormaps['viridis'])
	plt.contour(raster1, levels=15,lw=0.5, alpha=0.5)
	plt.scatter(coordsx,coordsy, s=150, color="black", facecolors='none', lw=3,alpha=1, label="Sites")
	plt.scatter(residencex1, residencey1, color="#7f0000", marker="+", lw=3, s=600, label="Residence")
	plt.scatter(residencex1, residencey1, color="red", marker="x", s=10)
	plt.scatter(residencex2, residencey2, color="#7f0000", marker="+", lw=3,s=600)
	plt.scatter(residencex2, residencey2, color="red", marker="x", s=10)
	
	
	global centerx, centery	
	plt.scatter(centerx, centery, color="blue", marker="x",label="Barycenter", s=600)
	
	drawmax=1
	
	if(drawmax==1):
		maxi1=np.max(raster1)
		indexes0=np.argwhere(raster1 == maxi1)
		indexes1=np.array(indexes0[0])
		print("I ", indexes1)
		print(indexes1[0])
		print(indexes1[1])	
		ix1=indexes1[1]
		iy1=indexes1[0]
		plt.scatter(ix1,iy1, color="#7f007f", marker="o",facecolors='none', lw=3, alpha=0.95, s=400, label="Maximum probability")
	
	
	
	
	plt.legend(fontsize=12)

	plt.xticks(fontsize=12)
	plt.yticks(fontsize=12)
	

	plt.tight_layout()
	plt.show()
	
	
	return(0)
	

def rossmo_1():
	global sizex, sizey, raster
	global coordsx,coordsy, numpoints
	global mindist
	global buf1, f, g,koo
	global standard_distance_manhattan
	global standard_distance_euclidian	
		
	disk=mindist/2
	print(disk)
	## official
	
	#buf1=disk*2
	#buf1=mindist/2
	buf1=standard_distance_manhattan/2
	#buf1=standard_distance_euclidian
	#buf1=buf1/sqrt(2)
	

	print("Rossmo analysis ")	
	print("center       ", centerx, centery)
	print("deviation    ",deviationx, deviationy)
	print("min distance ", mindist)
	print("buffer ",buf1)
	print("f      ", f)
	print("g      ", g);
	
	## tweaked
	#buf1=mindist
	## jn test
	
	
	#buf1=3
	
	##buf1=mindist
	
	
	
	for iy in range(0,sizey):
		for ix in range(0,sizex):
			pp=0
			for n in range(0,numpoints):
				px=coordsx[n]
				py=coordsy[n]
				dx=px-ix
				dy=py-iy
				dis=abs(dx)+abs(dy)
				difu=abs(dx)-abs(dy)
				
				fii=0
				if(dis==0): dis=1
				if(dis>buf1): fii=1
				###fii=0.5 ## testing exponent!
				term1=fii/(pow(dis,f))
				
				term2a1=(1-fii)*pow(buf1, (g-f))
				a2=(2*buf1)-abs(dx)-abs(dy)
				
				#if(a2==0): a2=1
				#print(a2)
				term2a2=pow( a2, g)
				#print(term2a1, term2a2)
				#if(term2a2<1): term2a2=1
				term2=term2a1/term2a2
				
				
				val1=koo*(term1+term2)
				
				#print(pp)
				
				p=val1
				
				pp=pp+p
							
			
			#print(pp)
			raster[iy,ix]=pp
			
					
		
	return(0)
	
def save_matrix_to_img():
	global raster
	nm1=normalize_matrix(raster)
	nm2=np.uint8(nm1*255)
	im1 = Image.fromarray(nm2,'L')
	#im1.show()
	im1.save("analysis1.png")
	return(0)
	

## main progg



## take account only n inner locations

#crime_locations_dropper(6) ## in lon lat coords mode ...
	
#quit(-1)
	

## if you use geographical coordinates
## translate them to raster coords
#translate_geographical_coordinates()




numpoints=len(coordsx)
	

mindist=find_nearest()


rossmo_1() ## ok	

#density_matrix_1()
##distance_matrix() ## mybe nok





#exponential_ring_1() ## ok 

#gaussian_ring_1() ##ok	

#gaussian_distance_from_center_1()
	

#exponential_distance_1()
	


## not ok

#gaussian_rossmooth_1() ## nok

## location of next probably attack
#gaussian_boundary_method() ## maybe OK location of next attack



#print(mindist)


save_matrix_to_img()

draw_raster_1()

## almost partially nok
save_netcdf("analysis1.nc")    
#draw_folium_map(raster) ## nok
Licensing

I, the copyright holder of this work, hereby publish it under the following license:
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
to share – to copy, distribute and transmit the work
to remix – to adapt the work
Under the following conditions:
attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history

Click on a date/time to view the file as it appeared at that time.
	Date/Time	Dimensions	User	Comment
current	10:50, 19 September 2023	1,047 × 757 (169 KB)	Merikanto (talk \| contribs)	Update of layout
	09:09, 18 September 2023	640 × 480 (118 KB)	Merikanto (talk \| contribs)	Update
	08:26, 15 September 2023	751 × 666 (172 KB)	Merikanto (talk \| contribs)	Update
	18:22, 14 September 2023	619 × 569 (146 KB)	Merikanto (talk \| contribs)	Uploaded own work with UploadWizard
You cannot overwrite this file.
File usage on Commons

There are no pages that use this file.
Metadata

This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Software used	Matplotlib version3.7.0, https://matplotlib.org/
Horizontal resolution	39.37 dpc
Vertical resolution	39.37 dpc
File:Rossmo analysis example desalvo murders 1.png

Captions

Captions

Summary

Licensing

File history

File usage on Commons

Metadata

Structured data

Items portrayed in this file

depicts

creator

some value

copyright status

copyrighted

copyright license

Creative Commons Attribution-ShareAlike 4.0 International

inception

14 September 2023

media type

image/png

source of file

original creation by uploader

Navigation menu

File:Rossmo analysis example desalvo murders 1.png

Captions

Captions

Summary

Licensing

File history

File usage on Commons

Metadata

Navigation menu

Search
