File:Rossmo analysis example desalvo murders 1.png
![File:Rossmo analysis example desalvo murders 1.png](https://upload.wikimedia.org/wikipedia/commons/thumb/1/19/Rossmo_analysis_example_desalvo_murders_1.png/800px-Rossmo_analysis_example_desalvo_murders_1.png?20230919105009)
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)
File information
Structured data
Captions
Captions
Rossmo analysis example: DeSalvo murders
Summary
editDescriptionRossmo analysis example desalvo murders 1.png |
English: Gaographic analysis of serial crime: attempt to guess location of residence of killer. |
Date | |
Source | Own work |
Author | Merikanto |
Sources of equations and data
-
- 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
-
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
editI, the copyright holder of this work, hereby publish it under the following license:
![w:en:Creative Commons](https://upload.wikimedia.org/wikipedia/commons/thumb/7/79/CC_some_rights_reserved.svg/90px-CC_some_rights_reserved.svg.png)
![attribution](https://upload.wikimedia.org/wikipedia/commons/thumb/1/11/Cc-by_new_white.svg/24px-Cc-by_new_white.svg.png)
![share alike](https://upload.wikimedia.org/wikipedia/commons/thumb/d/df/Cc-sa_white.svg/24px-Cc-sa_white.svg.png)
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 | Thumbnail | 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 | |
---|---|
Horizontal resolution | 39.37 dpc |
Vertical resolution | 39.37 dpc |