File:Atlanta murders 1979 1981 euclidian distance based analysis 1.png
![File:Atlanta murders 1979 1981 euclidian distance based analysis 1.png](https://upload.wikimedia.org/wikipedia/commons/thumb/8/81/Atlanta_murders_1979_1981_euclidian_distance_based_analysis_1.png/800px-Atlanta_murders_1979_1981_euclidian_distance_based_analysis_1.png?20230916084945)
Original file (1,280 × 816 pixels, file size: 1.42 MB, MIME type: image/png)
Captions
Captions
Summary
editDescriptionAtlanta murders 1979 1981 euclidian distance based analysis 1.png |
English: Atlanta murders 1979 1981 - euclidian distance based analysis |
Date | |
Source | Own work |
Author | Merikanto |
Used OpenStreetMap OSM
https://www.openstreetmap.org/copyright
CC BY-SA 2.0
Code
-
- 16.9.2023 v. 0000.0000.0021
-
-
- attempt to make geographical rossmo analysis of crime locations
- tries to shrikk area to find suspect of crime series
-
-
- atlanta 1979 or jack the ripper
-
- based on code at
-
-
-
-
- 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
- settings !
- rosmo matrix dims
dimx=100
dimy=100
- map initial zoom factor
zoomfac1=16
- draw
draw_suspect=1
- relative drawing skala of symbols
draw_skala=1
- extent of map
deltarectlon=0.12
deltarectlat=0.12
- symbols drawing scala
drawskala=1
- draw rossmo or another matrix
drawmatrix=1
- draw user image
drawimage=0
imagename1="ripper1.png"
- nw
- imglonnw1=-0.08173
- imglatnw1=51.52137
- imglonsw1=-0.08202
- imglatsw1=51.51208
- imglonne1=-0.05690
- imglatne1=51.52137
- imglonse1=-0.05679705085
- imglatse1=51.51195544111344
- imglonsw1=-0.08202
- imglatsw1=51.51208
- imglonne1=-0.05690
- imglatne1=51.52137
imglonsw1=0
imglatsw1=1
imglonne1=0
imglatne1=1
- crime locations, Atlanta
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]
- jack the ripper
- 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]])
- williams = [33.754065, -84.446577]
- names=["A","B","C","D","E"]
- lats<-c(51.5203785, 51.5186270, 51.5199207, 51.5137500,51.5137633)
- lons<-c(-0.0737000,-0.0752304,-0.0610239,-0.0779247,-0.0657004)
- dates<-as.Date(c('1888-01-01', '1888-01-01','1888-01-01','1888-01-01','1888-01-01') )
- lats=[51.5203785, 51.5186270, 51.5199207, 51.5137500,51.5137633]
- lons=[-0.0737000,-0.0752304,-0.0610239,-0.0779247,-0.0657004]
- names=["Nichols","Eddowes","Kelly","Chapman","Stride"]
- ages=[43,46,35,47,44]
- jack the ripper
- crimelocations=[[51.5203785,-0.0737000],[51.5186270,-0.0752304],[ 51.5199207, -0.0610239],[ 51.5137500,-0.0779247],[ 51.5137633, -0.0657004]]
- suspectlat = 51.5203785
- suspectlon = -0.0737000
- suspect= [suspectlat, suspectlon]
suspectx = 0
suspecty = 0
- 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)
- 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)
- second rossmo
- 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)
- 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),
html='
%s
' % "Max prob",
)
).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),
html='
%s
' % "Williams",
)
).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),
# html='
%s
' % "Barycenter",
# )
#).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),
# html='
%s
' % "Median",
# )
#).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),
# html='
%s
' % "MinMax",
# )
#).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)
- 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(5)
#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)
- create estimation of propability of location of suspect
calculate_probability(crimelocations, suspect)
Licensing
edit![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)
- 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 | 08:49, 16 September 2023 | ![]() | 1,280 × 816 (1.42 MB) | Merikanto (talk | contribs) | Update |
14:46, 4 September 2023 | ![]() | 1,280 × 816 (1.45 MB) | Merikanto (talk | contribs) | Update | |
17:44, 2 September 2023 | ![]() | 1,280 × 816 (990 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.