File:Kyllikki saari murder site rossmo 2 1 1 1.png

1. 1. https://github.com/hophopkhangaroo/rossmo

2. 1. rossmo formula with folium
3. 1. python 3.10 code anaconda windows
4. 30.8.2023 0000.0000.0000
5. 1. rossmo code based on
6. 1. https://github.com/hophopkhangaroo/rossmo

1. import point

import json import os import time import math import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sb from netCDF4 import Dataset,num2date,date2num from scipy.interpolate import griddata import scipy as sp import scipy.ndimage

from PIL import Image

import folium from folium import plugins from selenium import webdriver

import branca

import geojsoncontour

import geopandas as gpd

buffer = 150 f = 2/3 ## decay raye beyound buffer g = 3/4 ## decay rate within buffer

ncoutname="myfile.nc" varname="z"

1. 1. lat, lon rectangle

lon1=21.90 lon2=22.05 lat1=62.10 lat2=62.18

xsize=256*4 ysize=256*4

kenterlon1=21.951671 kenterlat1=62.144996

crime_geolocations=[ [21.951702795442856,62.14509129669217 ], [ 21.956429196411612,62.13349402990378 ] ]

sw1=[ lat1, lon1] ne1=[ lat2, lon2]

class Point():

   Creates coordinate points with values (x,y)

   def __init__(self, x, y):
       self.x = x
       self.y = y

   def taxicab_distance(self, other):
       Returns distance using taxicab metric
       dx = abs(self.x-other.x)
       dy = abs(self.y-other.y)
       return dx + dy

1. 1. WARNING not georeferenced netcdf output !!!

def ncout(ncoutname,varname, dataout1, lon1,lon2,lat1,lat2): #ny, nx = (ysize, xsize) ny, nx=np.shape(dataout1) lons = np.linspace(lon1,lon2,nx) lats = np.linspace(lat1, lat2,ny) #dataout = np.random.random((ny,nx)) ncout = Dataset(ncoutname,'w','NETCDF4') ncout.createDimension('lon',nx) ncout.createDimension('lat',ny) lonvar = ncout.createVariable('lon','float32',('lon')) #lonvar.long_name = 'lon' lonvar.units = 'degrees_east' lonvar.axis = 'lon' lonvar[:] = lons latvar = ncout.createVariable('lat','float32',('lat')) #latvar.long_name = 'lat' latvar.units = 'degrees_north' latvar.axis = 'lat' latvar[:] = lats myvar = ncout.createVariable(varname,'float32',('lon','lat')) #myvar.setncattr('units','degrees') #myvar.setncattr('grid_mapping', 'spatial_ref') #crs = ncout.createVariable('spatial_ref', 'i4') #crs.spatial_ref='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]' #crs.spatial_ref='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","4326"]]' myvar[:] = dataout1 ncout.close() return(0)

1. Grid Dimensions

nRows = xsize nCols = ysize

1. Known locations of committed crimes
2. crime_coords = [
3. Point(150*4,200*4),
4. Point(110*4,78*4),
5. ]

1. crime_locations = [
2. Point(392,510),
3. Point(390,428)
4. ]

crime_locations=[]

def calculate_radius_from_points(): ## experimental only auto-radius ... lenu=len(crime_locations) x00=crime_locations[0].x y00=crime_locations[0].y xr=0 yr=0 for n in range (1, lenu): x=crime_locations[n].x-x00 y=crime_locations[n].y-y00 xr=xr+x yr=yr+y

xr=xr/lenu yr=yr/lenu rr=int(math.sqrt(xr*xr+yr*yr)) print("Radius",rr) return(rr)

def view_map(rosmox):

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

imr=imr*2 imb=255-imb img=img*2.0 ima=255-imr ima2=ima/256

rosmoy=np.dstack((imr, img, imb, ima))

dimlat=lat2-lat1 dimlon=lon2-lon1

deltalon=dimlon/xsize deltalat=dimlat/ysize

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

for iy in range(0, ysize-1): tlat=lat2-deltalat*iy for ix in range(0, xsize-1): tlon=lon1+deltalon*ix tros=rosmox[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), xsize) lat_arr = np.linspace(np.min(lattab), np.max(lattab), ysize) 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=32 colors = ['black','black']

#levels = len(colors) cm = branca.colormap.LinearColormap(colors, vmin=vmin, vmax=vmax).to_step(levels) #contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.75, colors=colors, linestyles='None', vmin=vmin, vmax=vmax) contourf = plt.contourf(lon_mesh, lat_mesh, ros_mesh, levels, alpha=0.5, vmin=vmin, vmax=vmax) #contour = plt.contour(lon_mesh, lat_mesh, ros_mesh, alpha=0.2, colors=['black'], linewidths=[0.15]) geojson = geojsoncontour.contourf_to_geojson( contourf=contourf, min_angle_deg=0.01, ndigits=5, stroke_width=0.5, fill_opacity=0.2)

#kenterlon1=21.951671 #kenterlat1=62.144996 map = folium.Map(location=[kenterlat1, kenterlon1], TileProvider='OpenStreetMap', zoom_start=15, control_scale = True,) #folium.raster_layers.ImageOverlay(rosmoy,[[ lat1, lon1],[ lat2, lon2]],opacity=0.5,).add_to(map) #folium.LayerControl().add_to(map) # NOK #folium.GeoJson(geojson,).add_to(map) # # Plot the contour on Folium map #NOK folium.GeoJson( geojson, style_function=lambda x: { 'color': "#00007f", 'weight': 0.2, 'fillColor': None, 'opacity': 0.2, } ).add_to(map) #gdf = gpd.GeoDataFrame.from_features(json.loads(geojsoncontour.contour_to_geojson(contour=contour,min_angle_deg=3.0,ndigits=5,stroke_width=0.1))).set_crs("EPSG:4326") #map = gdf.explore(color=gdf["stroke"]) folium.raster_layers.ImageOverlay(rosmoy,[[ lat1, lon1],[ lat2, lon2]],opacity=0.5,).add_to(map) folium.LayerControl().add_to(map) # # # # latx1=crime_geolocations[0][1] lonx1=crime_geolocations[0][0] folium.Circle([latx1, lonx1], 45, fill_color="grey", opacity=1.0,color = 'grey', fill=True).add_child(folium.Popup('Hauta')).add_to(map)

latx1=crime_geolocations[1][1] lonx1=crime_geolocations[1][0] folium.Circle([latx1, lonx1], 45, fill_color="grey", opacity=1.0,color = 'grey', fill=True).add_child(folium.Popup('Hauta')).add_to(map)

map.fit_bounds([sw1, ne1]) delay=5 fn="testmap.html" tmpurl='file://{path}/{mapfile}'.format(path=os.getcwd(),mapfile=fn) map.save(fn) browser = webdriver.Firefox() browser.get(tmpurl) time.sleep(delay) browser.save_screenshot('testmap.png') browser.quit() return(0)

def print_crime_coords(): lenu=len(crime_locations) for n in (0, (lenu-1)): print(crime_locations[n].x,crime_locations[n].y ) return(0)

def add_crime_coordinates(): lenu=len(crime_geolocations) ulon=lon2-lon1 ulat=lat2-lat1 dlon=ulon/xsize dlat=ulat/ysize

for n in (0, (lenu-1)): lo=crime_geolocations[n][0] la=crime_geolocations[n][1] elon=lo-lon1 elat=la-lat1 dx=int(elon/dlon) dy=int(elat/dlat) crime_locations.append(Point(dx,dy))

return(0)

def main(): ## geo coords to relative grid coords ...

   add_crime_coordinates()
   buffer=calculate_radius_from_points()
   print("Buffer, f, g ",buffer, f , g)
   print_crime_coords()
   #quit(-1)
   aList = create_coord_list(nRows, nCols)
   aList = rossmo_list(aList)
   aTable = create_coord_table(aList)
   fig, ax = plt.subplots(figsize=(11, 9))
   #sb.heatmap(aTable, cmap="icefire")
   table1=np.array(aTable).astype(float)
   print (np.shape(table1))
   table2=np.flipud(table1)
   ncout("rosmos.nc","z", table2, lon1,lon2,lat1,lat2)
   #quit(-1)
   #plt.imshow(table1)
   #sb.heatmap(aTable, cmap="inferno")
   #ax.invert_yaxis()
   #plt.show()
   view_map(table2)
   return(0)
   

1. Creates a list to store cell coordinate locations

def create_coord_list(nRows, nCols):

   aList = list()
   for y in range(nRows):
       for x in range(nCols):
           aList.append(Point(x,y))
   return aList

1. Turns a list of cell coordinates into a table

def create_coord_table(aList):

   table = np.array(aList)
   table = table.reshape(nRows,nCols)
   table = pd.DataFrame(table)
   return table

1. Given a list of crime locations and a cell location, returns an array of all
2. distances from the crime

def crime_distance(cell, crime_locations):

   crime_distances = list()
   idx = 0
   for crime in crime_locations:
       crime_distances.append(cell.taxicab_distance(crime_locations[idx]))
       idx += 1
   return crime_distances

1. Returns an individual summand of Rossmo's equation

def rossmo_summand(distance, buffer, f, g):

   if distance > buffer:
       return 1/distance**f
   else:
       return (buffer**(f-g))/(2*buffer-distance)**g

1. Returns a likelihood for a cell using Rossmo's equation by summing over each
2. crime scene location

def rossmo_prob(crime_distances):

   px = 0
   for i in range(len(crime_distances)):
       if rossmo_summand(crime_distances[i], buffer, f, g) == None:
           px += 0
       else:
           px += rossmo_summand(crime_distances[i], buffer, f, g)
   return px

1. Given a list of cell coordinates, calculates likelihood for each cell

def rossmo_list(coords):

   aList = list()
   for coord in coords:
       crime_distances = crime_distance(coord, crime_locations)
       aList.append(rossmo_prob(crime_distances))
   return aList

main()

quit(-1)