* Redistributions of source code must retain the above copyright
     notice, this list of conditions and the following disclaimer.
   * Redistributions in binary form must reproduce the above copyright
     notice, this list of conditions and the following disclaimer in the
     documentation and/or other materials provided with the distribution.
   * Neither the name of Stamen Design nor the names of its contributors may
     be used to endorse or promote products derived from this software without
     specific prior written permission.

1. 1. exponent ring from crime locations
2. 1. R 4.30 script
3. 1. note slow in big raster sizes
4. 1. 15.9.2023 0000.0001bee

library(raster) library(sp)

1. 1. map view

library(leaflet) library(pandoc) library(htmlwidgets) library(webshot) library(mapview)

xsize=120 ysize=120

1. buffer1=xsize/20 ## radius of buffer

buffer1=0 ## note this resizes automati

f=0.5 ## inner exp coeff g=0.66 ## outer exp coeff

1. 1. nokf=0.5 ## inner exp coeff g=0.5 ## ouref exp coeff

1. f=0.33 g=0.66 ## nok

1. f=1 g=1 nok

crimelats<-c(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)

crimelons<- c(

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

1. 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]])

suspectlat <- c(33.754065) suspectlon <- c(-84.446577)

p1=c(-84.532406,33.703093) p2=c(-84.49509,33.660032) p3=c(-84.383959,33.741141) p4= c(-84.447227,33.711061) p5= c(-84.584169,33.701493) p6= c(-84.350482,33.746652) p7= c(-84.4941276,33.6605585) p8= c(-84.465294,33.755227) p9= c(-84.350066,33.692281) p10= c(-84.408613,33.738875) p11= c(-84.470401,33.805397) p12= c(-84.427292,33.677783) p13= c(-84.358048,33.679033) p14= c(-84.455166,33.858629) p15= c(-84.573247,33.68205) p16= c(-84.067554,33.68164) p17= c(-84.332364,33.837342) p18= c(-84.249387,33.680747) p19= c(-84.128966,33.631259) p20= c(-84.64159,33.683782) p21= c(-84.681008,33.653495) p22= c(-84.422132,33.766465) p23= c(-84.6803,33.653852) p24= c(-84.500141,33.802901) p25= c(-84.3287373,33.7392226) p26= c(-84.499154,33.804113)

centerlat1=mean(crimelats) centerlon1=mean(crimelons)

delta1=0.01

minlat1=min(crimelats)-delta1 maxlat1=max(crimelats)+delta1 minlon1=min(crimelons)-delta1 maxlon1=max(crimelons)+delta1

ext1<-c(minlon1, maxlon1, minlat1, maxlat1)

crs1 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

numpoints=length(crimelats)

pix1<-c() piy1<-c()

dimlon1=maxlon1-minlon1 dimlat1=maxlat1-minlat1

dellon1=dimlon1/xsize dellat1=dimlat1/ysize

for (n in 1:numpoints) { lona1=crimelons[n] lata1=crimelats[n] slona1=lona1-minlon1 slata1=lata1-minlat1

paax1=as.integer(slona1/dellon1) paay1=as.integer(slata1/dellat1)

pix1=append(pix1, paax1) piy1=append(piy1, paay1) }

1. print (pix1)
2. print(piy1)

1. 1. calculate minimum distance, if over 1 unit in raster

mindistance1=9999999

meandistancex=0 meandistancey=0; centerx=0 centery=0

for (m in 1:(numpoints))
{

px2=pix1[m] py2=piy1[m] centerx=centerx+px2 centery=centery+py2

for (n in 1:(numpoints))
{

px1=pix1[n] py1=piy1[n]

dx1=px1-px2 dy1=py1-py2

dd1<-sqrt(dx1*dx1+dy1*dy1)

if(dd1>1){ if(dd1<mindistance1) mindistance1=dd1 meandistancex=meandistancex+abs(dx1) meandistancey=meandistancey+abs(dy1)

}

}

}

centerx=centerx/numpoints centery=centery/numpoints

meandistancex=meandistancex/(numpoints*numpoints-numpoints) meandistancey=meandistancey/(numpoints*numpoints-numpoints)

meandistance1=sqrt(meandistancex*meandistancex+meandistancey*meandistancey)

print(numpoints) print (centerx) print (centery) print(mindistance1) print(meandistance1)

1. buffer1=mindistance1

1. buffer1=mindistance1*sqrt(2) ## ok OK

1. 1. buffer1=mindistance1*2 ## nok

1. buffer1=meandistance1/sqrt(2) ## quite ok

buffer1=sqrt(meandistance1) ## ok

1. buffer1=meandistance1/2

print("Buffer:") print(buffer1)

1. stop(-1)

ras1<-raster(nrows=xsize, ncols=ysize)

values(ras1) <- runif(ncell(ras1))*0

extent(ras1) <- c(0,xsize,0,ysize)

1. numpoints=5

1. pix1<-runif(numpoints)*xsize
2. piy1<-runif(numpoints)*ysize

for (iy in 1:(ysize)) { for (ix in 1:(xsize)) {

valu1=0


for (n in 1:(numpoints))
{

tx=pix1[n] ty=piy1[n]

#pax=ras1[iy,ix] dx=abs(tx-ix) dy=abs(ty-iy) #print (iy) #print(dx) dist1=sqrt(dx*dx+dy*dy) ## eucleidian #dist1=dx+dy ## manhattan

1. if(dist1>buffer1)
2. {
3. a1=exp(-dist1/buffer1)
4. }

dipu1=dist1-buffer1 adipu1=abs(dipu1) redipu1=dipu1/buffer1 ## relative radius vs buffer aredipu1=adipu1/buffer1 ##relative abs radius vs buffer ## normal #a1=pnorm(redipu1, 1, 0.5)

## exponent #a1=exp(-redipu1) if(redipu1>1) a1=exp(-aredipu1/g) if(redipu1<=1) a1=exp(-aredipu1/f)

## 1:st point is first, wthis is weighted point ##

#w=3 ## weight of points comes from order #w=10/sqrt(n) ##w=1 ##if(n==1) w=3 w=1 ## all points equal

a2=a1*w

#print(dist1) valu1=valu1+a2 }

ras1[iy,ix]=valu1 

}

1. stop(-1)

1. print(iy)

}

min1<-minValue(ras1) max1<-maxValue(ras1) del1<-max1-min1

ras2<-((ras1-min1)/del1) rout1<-flip(ras2*100)

1. plot(flip(ras2))

1. points(pix1,piy1)

crs(rout1)<-crs1 extent(rout1)<-ext1

plot(flip(rout1))

points(crimelons, crimelats)

points(suspectlon, suspectlat, pch = 12, cex=2, col="black", bg="red", lwd=2)

f <- writeRaster(rout1, "estimation1.tif", format="Gtiff", overwrite=TRUE)

idx1 = which.max(rout1) pos1 = xyFromCell(rout1,idx1)

print("max pos") print(pos1)

amaaxlon1=as.numeric(pos1[1]) amaaxlat1=as.numeric(pos1[2])

print (amaaxlon1) print (amaaxlat1)

1. stop(-1)

1. Stadia.StamenTonerBackground

1. tiles1="https://stamen-tiles.a.ssl.fastly.net/toner/{z}/{x}/{y}.png"

tiles1="https://tile.openstreetmap.org/{z}/{x}/{y}.png"

1. tiles1="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png" ## satellite
2. tiles1="https://{s}.tile-cyclosm.openstreetmap.fr/cyclosm/{z}/{x}/{y}.png"
3. tiles1="https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}.png"
4. tiles1="https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}.png"
5. tiles1="https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png"
6. tiles1="https://{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}{r}.png"
7. tiles1="https://{s}.basemaps.cartocdn.com/dark_nolabels/{z}/{x}/{y}{r}.png"
8. tiles1=https://{s}.basemaps.cartocdn.com/rastertiles/voyager_nolabels/{z}/{x}/{y}{r}.png"

1. tiles1="https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"

tiles1="https://tiles.stadiamaps.com/tiles/stamen_toner_background/{z}/{x}/{y}{r}.png" #

1. tiles1="https://tiles.stadiamaps.com/tiles/stamen_toner_lines/{z}/{x}/{y}{r}.png"

contours1 <- rasterToContour(rout1)

pal = colorNumeric(c("blue", "green", "yellow", "orange", "red", "violet", "violet"), 50:100,na.color = "transparent")

map1 <- leaflet(padding=0) %>%

 addTiles(tiles1) %>% 
 fitBounds(minlon1,minlat1,maxlon1,maxlat1) %>% 
 setView(centerlon1, centerlat1,  zoom = 12) %>%
 addRasterImage(rout1, colors = pal, opacity = 0.5) %>%
 #addRasterImage(rout1, opacity = 0.5) %>%
 addPolylines(data =contours1,fillOpacity=0.5,fillColor = "transparent",opacity=0.5,weight=1) %>%
 addCircleMarkers(lng = crimelons, lat =  crimelats, radius = 8, color = "#F00", popup = NULL) %>%
 fitBounds(minlon1,minlat1,maxlon1,maxlat1) %>%
 addCircleMarkers(suspectlon, suspectlat, radius = 10, color = "#003", popup = NULL) %>%
 addCircleMarkers(centerlon1, centerlat1, radius = 10, color = "#0f0", popup = NULL)  %>%
 addCircleMarkers(amaaxlon1, amaaxlat1, radius = 10, color = "#00f", popup = NULL)  

1. addMarkers(suspectlon, suspectlat)
2. addMarkers(crimelons, crimelats)

1. map1<-fitBounds(minlon1,minlat1,maxlon1,maxlat1)

saveWidget(map1, file="map1.html", selfcontained = FALSE)

1. webshot("map1.html", file = "map1.png",cliprect = "viewport", remove_controls = c("zoomControl", "layersControl", "homeButton", "scaleBar",
2. "drawToolbar", "easyButton") )

1. mapshot(map1, file = paste0(getwd(), "/map1.png"))

willremove1<-c("zoomControl", "layersControl", "homeButton", "scaleBar","drawToolbar", "easyButton")

mapshot(map1, file = paste0(getwd(), "./map1.png"),remove_controls = willremove1)

print(".")