File:Exponential ring geographical analysis of locations sample 1.png

1. include <iostream>
2. include <fstream>
3. include <cstdint>

//#include <mem.h> // Linux

1. include <string.h> // Windows
2. include <math.h>

1. include "libattopng.h"

// exponent ring params

float buffer1=100;

//float f=0.500; //float g=0.667;

float f=1.500; float g=1.50;

int max_val=255; int sizex=1000; int sizey=1000;

// bounding box dimsonson from its center degrees double delta1=0.5;

int centerx, centery=0;

double mindistance=99999999; double meandistance=0;

int numpoints=26;

int pointsx[256]={600,400,800}; int pointsy[256]={300,700, 900};

double crimelats[256]={ 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 };

double crimelons[256] ={

-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 };

double suspectlat = 33.754065 ; double suspectlon = -84.446577 ;

double minlat1=suspectlat-delta1; double maxlat1=suspectlat+delta1; double minlon1=suspectlon-delta1; double maxlon1=suspectlon+delta1;

unsigned char *image1=NULL;

int setvalues(void) {

   int n=0;
   int x,y=0;
   double klat1,klon1, dlat1, dlon1=0;
   double px,py,qx,qy=0;
   double dx,dy,cx, cy,ex,ey=0;
   double diss1;

   minlat1=suspectlat-delta1;
   maxlat1=suspectlat+delta1;
   minlon1=suspectlon-delta1;
   maxlon1=suspectlon+delta1;

   printf("\nBounding box rectangle of output raster: ");
   printf("\nlon lat (%f %f) - (%f %f)",minlon1, minlat1, maxlon1, maxlat1);
   printf("\n\n R-style extent of output raster: ");
   printf("\n extent1<-c(%f,%f,%f,%f)",minlon1,maxlon1, minlat1,  maxlat1);

   double latcoeff1=sizey/(maxlat1-minlat1);

   double loncoeff1=sizex/(maxlon1-minlon1);

   printf("\n\n Points lon, lat ---> x, y\n");
   for (n=0;n<numpoints;n++)
       {
       klon1=crimelons[n];
       klat1=crimelats[n];

       dlat1=klat1-minlat1;
       dlon1=klon1-minlon1;

       x=int(dlon1*loncoeff1);
       y=int(dlat1*latcoeff1);

       printf("\n %f,%f ---> %i,%i",klon1, klat1, x, y);

       //printf("\n %f %f  %i %i",klon1, klat1, x, y);

       pointsx[n]=x;
       pointsy[n]=y;
   }

   printf("\n\nPython coordinates of input locations in raster units:");

printf("\ncoordsx=np.array([");

    for (n=0;n<numpoints;n++)
       {
       if(n<(numpoints-1)) printf("%i,",pointsx[n]) ;
       else printf("%i",pointsx[n]) ;

       }

printf("\n])");

 printf("\ncoordsy=np.array([");

    for (n=0;n<numpoints;n++)
       {
       if(n<(numpoints-1)) printf("%i,",pointsy[n]) ;
       else printf("%i",pointsy[n]) ;
       }

printf("\n])");

cx=0; cy=0;

int m=0;

   for (m=0;m<numpoints;m++)
    {

qx=pointsx[m]; qy=pointsy[m]; cx=cx+qx; cy=cy+qy;

for (n=0;n<numpoints;n++) { if(n!=m) { px=pointsx[n]; py=pointsy[n];

dx=qx-px; dy=qy-py;

ex=ex+fabs(dx); ey=ey+fabs(dy);

diss1=sqrt(dx*dx+dy*dy);

if(diss1<mindistance) { if(diss1>1) mindistance=diss1; }

}

} }

centerx=int(cx/numpoints);
centery=int(cy/numpoints);

ex=ex/(numpoints*numpoints-numpoints);
ey=ey/(numpoints*numpoints-numpoints);

meandistance=sqrt(ex*ex+ey*ey);

printf("\nmean distance %f", meandistance);

//buffer1=meandistance/sqrt(2);

//   buffer1=sqrt(meandistance);

buffer1=meandistance/2;

printf("\n Buffer %f", buffer1);

return(0); }

1. define RGBA(r, g, b, a) ((r) | ((g) << 8) | ((b) << 16) | ((a) << 24))

int savepng(char *filename1, unsigned char *data1, int rows, int cols, int max_val) { int x, y,gee;

//libattopng_t* png = libattopng_new(sizex, sizey, PNG_RGBA); libattopng_t* png = libattopng_new(cols, rows, PNG_GRAYSCALE);

/*

   png = libattopng_new(W, H, PNG_GRAYSCALE);

   libattopng_start_stream(png, 0, 0);
   for (x = 0; x < W * H; x++) {
       libattopng_put_pixel(png, 255 * x / (W * H));
   }
   libattopng_save(png, "test_gray_stream.png");
   libattopng_destroy(png);

• /

 libattopng_start_stream(png, 0, 0);

for (y = 0; y < sizey; y++) {

 for (x = 0; x < sizex; x++) {
       gee=data1[y*cols+x];

// libattopng_set_pixel(png, x, y, RGBA(gee & 255, gee & 255, gee & 255, 255 ) );
   libattopng_put_pixel(png, gee);
 }

} libattopng_save(png, filename1); libattopng_destroy(png); return(0); }

int savepgm(char *filename1, unsigned char *data1, int rows, int cols, int max_val) {

   int i, j,pixel=0;
   int pixloc=0;
   std::fstream newfile (filename1, std::ios::out| std::ios::binary | std::ios::ate );

  // newfile.open(filename1, std::ios::app | std::ios::binary );
  // newfile.open(filename1);

   newfile << "P5" << std::endl;
   newfile << cols << " " << rows << std::endl;
   ///newfile << max_val << std::endl;
   newfile << max_val;
   for (  j = 0; j < rows; j++ )
   {
       for ( i = 0; i < cols; i++ )
       {
           pixloc=j*cols+i;
           pixel = data1[pixloc];
           newfile << (unsigned char)pixel;
       }

   }

    newfile.close();

return(0);

}

int free_image(unsigned char *imagee) {

   if(imagee)
   {
       free(imagee);
       imagee=NULL;
   }
return(0);

}

unsigned char * generate_image(int rows, int cols) {

   double pp,p=0;

   image1=(unsigned char *)malloc(rows*cols*4);

   if(!image1) return(NULL);

   memset(image1, 0,rows*cols*4 );

    for ( int j = 0; j < rows; j++ )
   {
       for ( int i = 0; i < cols; i++ )
       {
       pp=0;
         for ( int m = 0; m < numpoints; m++ )
           {

           int x=pointsx[m];
           int y=pointsy[m];
           int dx=i-x;
           int dy=j-y;
     //      int adx=abs(dx);
       //    int ady=abs(dy);

           double dist1=double(sqrt(dx*dx+dy*dy));
        //   int mdist1=adx+ady;

         //  int idist1=0;

           if(dist1==0) dist1=1;

           float reldist1=(float)(dist1-buffer1)/(float)buffer1;

           float absreldist1=fabs(reldist1);

           float fun1=0;

           double funk1=256/numpoints;

           if(reldist1>1) fun1=funk1*exp(-g*absreldist1);
           if(reldist1<=1) fun1=funk1*exp(-f*absreldist1);

           p=fun1;

           pp=pp+p;

           }
               image1[j*cols+i]=int(pp);

       }
   }

return(image1);

}

int normalize(unsigned char*image) { int n, imalen1=0; int minima, maxima=00;

minima=255; maxima=0; imalen1=sizex*sizey;

for(n=0;n<imalen1;n++) {

}

}

int main (int argc, char **argv) {

   setvalues();

   image1=generate_image(sizey, sizex);

   savepgm((char *)"outp.pgm", image1,  sizey,  sizex, 255);
   savepng((char *)"out.png", image1,  sizey,  sizex, 255);

   free_image(image1);

}

2. 1. view external crime locations geographic exponential ring analysis raster
4. R 4.3.0, leaflet, stamen toner no symbols ...
5. 15.9.2023 0000.0001b

library(raster) library(rgdal)

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

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

iname1="out.png" oname1="analysis1.tif"

1. 1. atlanta serial killer 1979-1981

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)

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

minlon1=extent1[1] minlat1=extent1[3] maxlon1=extent1[2] maxlat1=extent1[4]

1. print(minlon1)

1. stop(-1)

1. 1. load external raster

rin1<-raster(iname1)

rout00<-flip(rin1)

max1=maxValue(rout00) min1=minValue(rout00) del1=max1-min1

rout01<-(rout00-min1)/del1

rout01=rout01*100

1. 1. small values out

threshold1=75

rout01=rout01-threshold1

rout01[rout01<0]=0

max1=maxValue(rout01) min1=minValue(rout01) del1=max1-min1

rout0<-(rout01-min1)/del1

rout0=rout0*100

crs(rout0) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

extent(rout0)<-extent1

f <- writeRaster(rout0, "analysis1.tif", format="GTiff", overwrite=TRUE)

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

print("max pos") print(pos1)

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

print (amaaxlon1) print (amaaxlat1)

1. 1. draw map1

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

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

contours1 <- rasterToContour(rout0, nlevels=10)

map1 <- leaflet() %>%

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

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

1. 1. save output map to disc

willremove1<-c("zoomControl", "layersControl", "homeButton", "scaleBar",

   "drawToolbar", "easyButton")

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

print(".")

1. 1. rossmo etc image visualization
2. 3 Python 3 code
4. 15.9.2023 0000.0000

import numpy as np import matplotlib.pyplot as plt from PIL import Image

sizey=1000 sizex=1000

1. lon

1. lat -84.946577 33.254065 -83.946577 34.254065

1. 1. 800x800
2. coordsx=np.array([331,361,450,399,289,476,361,385,477,430,380,415,470,393,298,703,491,557,654,243,212,419,213, 357, 494,357])
3. coordsy=np.array([359,324,389,365,357,394,325,400,350,387,441,338,339,483,342,342,466,341,301,343,319,409,319,439,388,440])

coordsx=np.array([414,451,562,499,362,596,452,481,596,537,476,519,588,491,373,879,614,697,817,304,265,524,266,446,617,447]) coordsy=np.array([449,405,487,456,447,492,406,501,438,484,551,423,424,604,427,427,583,426,377,429,399,512,399,548,485,550])

coordsy=sizey-coordsy

1. delta=0.3

1. pic0 = Image.open("./out.png").convert('RGB')

pic0 = Image.open("./out.png").convert('L')

arr1=np.flipud(np.array(pic0))

1. arr1=np.array(pic0)

1. Creates PIL image
2. img = Image.fromarray(pic0, 'L')

1. pic0.show()

min1=np.min(arr1) max1=np.max(arr1) del1=max1-min1

arr2=1*(arr1-min1)/del1

plt.imshow(arr2) c = plt.contour(arr2, levels=[ 0.9,0.95, 0.98], colors=["orange", "red", "violet"], alpha=0.5)

plt.scatter(500,500, color='red', s=32, marker="x") plt.scatter(coordsx, coordsy, color='blue', s=12)

drawmax=1

if(drawmax==1): maxi1=np.max(arr2) indexes0=np.argwhere(arr2 == 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=1, alpha=1, s=200, label="Maximum probability")

plt.show()