Search code examples
raggregationraster

Aggregation of coordinates in a grid


I'm using these coordinates (raster) in order to plot them above an empty image.

rm(list = setdiff(ls(), lsf.str()))
mydata=structure(list(Position_X = c("95", "80", "43", "29", "92", "106", 
"95", "105", "80", "78", "74", "109", "91", "125", "110", "46", 
"48", "57", "54", "34", "76", "136", "82", "75", "130", "132", 
"134", "136", "130", "130", "78", "93", "79", "127", "113", "113", 
"128", "117", "97", "80", "132", "88", "114", "51", "122", "97", 
"120", "85", "103", "154", "161", "165", "164", "207", "185", 
"207", "165", "193", "184", "160", "247", "207", "208", "209", 
"230", "211", "154", "168", "214", "189", "186", "207", "242", 
"197", "200", "163", "214", "189", "246", "217", "167", "199", 
"167", "186", "162", "161", "201", "186", "134", "221", "146", 
"177", "164", "158", "163", "249", "136", "176", "257", "308", 
"357", "372", "295", "311", "342", "309", "360", "355", "311", 
"312", "370", "316", "358", "1450", "352", "337", "378", "340", 
"296", "334", "342", "308", "312", "312", "303", "374", "267", 
"327", "340", "376", "389", "341", "368", "262", "353", "343", 
"338", "354", "311", "352", "358", "356", "302", "406", "448", 
"452", "476", "392", "422", "480", "405", "416", "458", "406", 
"487", "480", "437", "411", "450", "464", "449", "467", "477", 
"484", "406", "406", "487", "473", "464", "435", "418", "470", 
"480", "504", "515", "467", "441", "422", "422", "380", "474", 
"428", "481", "482", "415", "497", "479", "419", "405", "385", 
"449", "479", "488", "477", "479", "364", "412", "378", "572", 
"568", "555", "566", "557", "556", "522", "556", "513", "555", 
"558", "564", "558", "555", "541", "555", "534", "557", "557", 
"557", "558", "522", "525", "519", "520", "569", "556", "516", 
"581", "540", "512", "535", "536", "556", "523", "533", "541", 
"557", "523", "537", "557", "524", "542", "517", "558", "509", 
"543", "590", "618", "660", "706", "697", "618", "679", "646", 
"639", "670", "625", "660", "620", "697", "636", "667", "703", 
"686", "679", "730", "728", "661", "692", "674", "677", "675", 
"624", "624", "696", "661", "712", "626", "683", "695", "672", 
"621", "700", "672", "632", "610", "631", "683", "636", "647", 
"629", "655", "688", "614", "663", "795", "753", "811", "786", 
"793", "850", "841", "788", "822", "810", "794", "808", "734", 
"792", "815", "731", "821", "805", "807", "809", "837", "809", 
"835", "839", "835", "782", "841", "839", "850", "852", "741", 
"824", "830", "843", "808", "823", "824", "787", "872", "848", 
"839", "810", "813", "790", "819", "957", "919", "846", "919", 
"903", "924", "867", "843", "861", "953", "925", "869", "908", 
"899", "1016", "1060", "1075", "1084", "1035", "992", "977", 
"958", "975", "908", "913", "962", "1063", "1008", "908", "958", 
"1056", "1013", "860", "851", "1001", "909", "891", "1034", "953", 
"958", "916", "953", "1021", "1027", "1035", "906", "956", "980", 
"1017", "1061", "925", "1044", "1051", "1050", "1048", "918", 
"1042", "1052", "1053", "1063", "1030", "1051", "1049", "998", 
"959", "1048", "979", "1003", "984", "1002", "1048", "992", "1016", 
"1042", "1039", "874", "1033", "917", "958", "863", "1048", "1107", 
"1104", "1105", "1133", "1152", "1109", "1009", "953", "1120", 
"1178", "1154", "1102", "976", "1057", "1191", "1195", "955", 
"1022", "986", "958", "621", "864", "882", "1106", "1177", "1197", 
"1062", "898", "1033", "916", "1115", "1196", "1057", "1103", 
"1058", "1170", "1168", "921", "917", "1151", "1129", "1197", 
"1065", "1156", "1160", "1115", "1121", "1176", "1108", "1127", 
"1073", "1175", "1085", "1419", "1156", "1154", "1375", "1376", 
"1508", "1126", "1341", "1441", "1470", "1381", "1325", "1462", 
"1463", "1360", "1262", "1505", "1404", "1178", "1177", "1395", 
"1373", "1313", "1386", "1388", "1518", "1285", "1248", "1229", 
"1181", "1465", "1328", "1103", "1392", "1501", "1465", "1410", 
"1397", "1375", "1328", "1456", "1459", "1323", "1159", "1166", 
"1508", "1303", "1329", "1454", "1459", "1147", "1228", "1421", 
"1316", "1456", "1340", "1341", "1397", "1447", "1309", "1400", 
"1483", "1253", "1303", "1296", "1302", "1299", "1150", "1228", 
"1470", "1526"), Position_Y = c("46", "43", "50", "35", "75", 
"78", "38", "44", "64", "62", "65", "33", "34", "52", "44", "44", 
"44", "43", "64", "58", "63", "43", "63", "36", "70", "68", "68", 
"68", "68", "66", "63", "80", "67", "40", "80", "75", "48", "59", 
"50", "41", "43", "70", "74", "33", "73", "78", "77", "39", "73", 
"51", "42", "84", "48", "53", "72", "45", "77", "72", "87", "48", 
"52", "55", "54", "76", "75", "74", "47", "73", "50", "52", "79", 
"48", "80", "81", "51", "77", "51", "53", "81", "50", "52", "52", 
"74", "78", "44", "83", "44", "77", "82", "80", "50", "73", "53", 
"46", "81", "86", "50", "53", "54", "80", "67", "82", "68", "89", 
"83", "85", "59", "87", "56", "51", "94", "43", "86", "58", "85", 
"56", "98", "44", "55", "92", "93", "76", "86", "52", "85", "96", 
"87", "50", "55", "96", "63", "91", "51", "51", "90", "86", "87", 
"66", "74", "86", "60", "60", "61", "91", "98", "96", "53", "57", 
"97", "51", "56", "95", "95", "57", "52", "97", "53", "95", "95", 
"52", "96", "96", "96", "52", "57", "59", "52", "97", "96", "53", 
"95", "53", "96", "52", "103", "53", "102", "96", "52", "101", 
"97", "98", "96", "95", "43", "51", "96", "48", "57", "63", "47", 
"46", "53", "96", "46", "55", "57", "63", "87", "61", "81", "74", 
"71", "81", "74", "72", "95", "77", "85", "61", "81", "82", "66", 
"72", "59", "70", "81", "84", "88", "73", "73", "55", "73", "85", 
"62", "81", "56", "57", "82", "82", "83", "54", "83", "60", "82", 
"68", "87", "90", "66", "91", "82", "52", "71", "98", "65", "65", 
"57", "54", "49", "97", "92", "99", "55", "50", "97", "94", "96", 
"55", "98", "95", "55", "49", "48", "102", "96", "102", "55", 
"96", "55", "95", "56", "56", "97", "52", "56", "101", "56", 
"100", "52", "96", "94", "101", "102", "49", "91", "98", "101", 
"49", "55", "50", "50", "98", "92", "95", "86", "54", "83", "67", 
"86", "86", "65", "66", "81", "65", "88", "86", "54", "92", "91", 
"96", "94", "58", "65", "62", "99", "72", "80", "69", "93", "88", 
"83", "66", "87", "85", "56", "66", "70", "68", "88", "68", "69", 
"60", "63", "89", "63", "67", "63", "66", "62", "59", "90", "90", 
"70", "98", "85", "88", "96", "64", "84", "85", "63", "79", "59", 
"58", "72", "63", "57", "91", "53", "98", "85", "87", "81", "70", 
"86", "83", "96", "72", "70", "85", "97", "94", "68", "95", "94", 
"72", "95", "83", "67", "90", "76", "62", "98", "60", "82", "58", 
"63", "93", "68", "95", "55", "62", "63", "66", "86", "60", "85", 
"59", "70", "60", "65", "63", "103", "91", "104", "64", "96", 
"96", "94", "62", "106", "63", "63", "63", "88", "71", "61", 
"69", "88", "93", "67", "91", "82", "66", "86", "92", "60", "94", 
"99", "102", "96", "62", "92", "67", "100", "53", "62", "62", 
"85", "62", "60", "60", "95", "59", "61", "50", "93", "93", "59", 
"61", "80", "103", "63", "95", "91", "71", "98", "95", "80", 
"61", "68", "56", "95", "57", "65", "99", "57", "101", "59", 
"67", "62", "60", "60", "96", "96", "97", "51", "93", "54", "63", 
"97", "56", "67", "81", "78", "91", "62", "51", "88", "64", "59", 
"65", "94", "54", "50", "77", "49", "104", "94", "61", "63", 
"65", "64", "91", "68", "61", "93", "91", "88", "96", "50", "104", 
"68", "64", "86", "67", "63", "94", "65", "67", "55", "90", "56", 
"50", "104", "100", "68", "55", "50", "99", "52", "98", "78", 
"50", "60", "60", "60", "75", "78", "69", "94", "70", "100", 
"55")), .Names = c("Position_X", "Position_Y"), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 
16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 
29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 
42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 
55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 
68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 
81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 
94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L, 
106L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 
117L, 118L, 119L, 120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L, 
128L, 129L, 130L, 131L, 132L, 133L, 134L, 135L, 136L, 137L, 138L, 
139L, 140L, 141L, 142L, 143L, 144L, 145L, 146L, 147L, 148L, 149L, 
150L, 151L, 152L, 153L, 154L, 155L, 156L, 157L, 158L, 159L, 160L, 
161L, 162L, 163L, 164L, 165L, 166L, 167L, 168L, 169L, 170L, 171L, 
172L, 173L, 174L, 175L, 176L, 177L, 178L, 179L, 180L, 181L, 182L, 
183L, 184L, 185L, 186L, 187L, 188L, 189L, 190L, 191L, 192L, 193L, 
194L, 195L, 196L, 197L, 198L, 199L, 200L, 201L, 202L, 203L, 204L, 
205L, 206L, 207L, 208L, 209L, 210L, 211L, 212L, 213L, 214L, 215L, 
216L, 217L, 218L, 219L, 220L, 221L, 222L, 223L, 224L, 226L, 227L, 
228L, 229L, 230L, 231L, 232L, 233L, 234L, 235L, 236L, 237L, 238L, 
239L, 240L, 241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L, 
250L, 251L, 252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 260L, 
261L, 262L, 263L, 264L, 265L, 266L, 267L, 268L, 269L, 270L, 271L, 
272L, 273L, 274L, 275L, 276L, 277L, 278L, 279L, 280L, 281L, 282L, 
283L, 284L, 285L, 286L, 287L, 288L, 289L, 290L, 291L, 292L, 293L, 
294L, 295L, 296L, 297L, 298L, 299L, 300L, 301L, 302L, 303L, 304L, 
305L, 306L, 307L, 308L, 309L, 310L, 311L, 312L, 313L, 314L, 315L, 
316L, 317L, 318L, 319L, 320L, 321L, 322L, 323L, 324L, 325L, 326L, 
327L, 328L, 329L, 330L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 
338L, 339L, 340L, 341L, 342L, 343L, 344L, 345L, 346L, 347L, 348L, 
349L, 350L, 351L, 352L, 353L, 354L, 355L, 356L, 357L, 358L, 359L, 
360L, 361L, 362L, 363L, 364L, 365L, 366L, 367L, 368L, 369L, 370L, 
371L, 372L, 373L, 374L, 375L, 376L, 377L, 378L, 379L, 380L, 381L, 
382L, 383L, 384L, 385L, 386L, 387L, 388L, 389L, 390L, 391L, 392L, 
393L, 394L, 395L, 396L, 397L, 398L, 399L, 400L, 401L, 402L, 403L, 
404L, 405L, 406L, 407L, 408L, 409L, 410L, 411L, 412L, 413L, 414L, 
415L, 416L, 417L, 418L, 419L, 420L, 421L, 422L, 423L, 424L, 425L, 
426L, 427L, 428L, 429L, 430L, 431L, 432L, 433L, 434L, 435L, 436L, 
437L, 438L, 439L, 440L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 
448L, 449L, 450L, 451L, 452L, 453L, 454L, 455L, 456L, 457L, 458L, 
459L, 460L, 461L, 462L, 463L, 464L, 465L, 466L, 467L, 468L, 469L, 
470L, 471L, 472L, 473L, 474L, 475L, 476L, 477L, 478L, 479L, 480L, 
481L, 482L, 483L, 484L, 485L, 486L, 487L, 488L, 489L, 490L, 491L, 
492L, 493L, 494L, 495L, 496L, 497L, 498L, 499L, 500L, 501L, 502L, 
503L, 504L, 505L, 506L, 507L, 508L, 509L, 510L, 511L, 512L, 513L, 
514L, 515L, 516L, 517L, 518L, 519L, 520L, 521L, 522L, 523L, 524L, 
525L, 526L, 527L, 528L, 529L, 530L, 531L, 532L, 533L, 534L, 535L, 
536L, 537L, 538L, 539L, 540L, 541L, 542L, 543L, 544L, 545L), class = "data.frame")

library(ggplot2)
library(png)
library(dplyr)
setwd("/Users/Folder/")
resize.win <- function(Width=6, Height=6)
{
    dev.off();
    windows(record=TRUE, width=Width, height=Height)
}
resize.win(1552,300)
#Replace the directory and file information with your info
img <- readPNG("plan3.png")
height<-dim(img)[1]
width<-dim(img)[2]
mypath <- file.path("C:/Users/Folder/img)
png(filename = mypath, width = 1552, height = 300)
#Set up the plot area
par(bty = 'n') 
plot(1:1, type='n', main="", xlab="", ylab="", xaxt='n', yaxt='n', ann=FALSE, xlim=c(1,1552),ylim=c(1,112))
#recuperation des infos
lim <- par()
rasterImage(img,0,0,width,height)
lines(mydata$Position_X, mydata$Position_Y, col = "red" ,type="p", pch=16,lwd=2, cex = 1.8)
par(xpd=NA,oma=c(17.9,25,11,10)) 
legend(par("usr")[1],par("usr")[3],legend = "red" ,text.font=2, fill="red", ncol = 3,bty = "n", cex = 0.8,lty=0,xjust=0, yjust=1.8)
dev.off()

The result is here :

enter image description here

Now, I want to aggregate these coordinates in order to have a more macroscopic representation, let's say every 3px*3px there's aggregation of points (coordinate) with different colors of aggregation's number.

It is possible to do this, and how ? Thanks a lot.

EDIT 1 : Use of DBSCAN

 mydata$Position_X=as.numeric(mydata$Position_X)
    mydata$Position_Y=as.numeric(mydata$Position_Y)
    df=mydata
    library("fpc")
    # Compute DBSCAN using fpc package
    set.seed(123)
    db <- fpc::dbscan(df, eps = 0.15, MinPts = 5)
    # Plot DBSCAN results
    plot(db, df, main = "DBSCAN", frame = FALSE)

Problem : This code only make plot of my points, no aggregation. Why I couldn't use mydataset in this way ?

EDIT2 : Raster

    coordinates(mydata) = ~Position_X+Position_Y
library(raster)
r <- raster(xmn=0, ymn=0, xmx=416, ymx=73, res=5)
r[] <- 0
xy <- mydata
tab <- table(cellFromXY(r, xy))
r[as.numeric(names(tab))] <- tab
plot(r)
points(xy, pch=20)

enter image description here


Solution

  • I write directly what I had in mind with raster::rasterize, which does directly fill a raster with the counts. Note that I used a resolution of 10 pixels, but you can choose anything else.

    # Libraries ----
    library(sp)
    library(raster)
    library(dplyr)
    library(ggplot2)
    library(rasterVis)
    
    # Convert data to numeric ----
    mydata <- mydata %>%
      mutate(Position_X = as.numeric(Position_X),
             Position_Y = as.numeric(Position_Y))
    
    # ==== With raster ====
    mydata.sp <- mydata
    coordinates(mydata.sp) <- ~Position_X+Position_Y
    
    r <- raster(xmn = xmin(mydata.sp), ymn = ymin(mydata.sp),
                xmx = xmax(mydata.sp), ymx = ymax(mydata.sp),
                res = 10)
    r[] <- 0
    plot(r)
    points(mydata.sp, pch = 20)
    
    # raster with classical plot
    mydata.r <- rasterize(mydata.sp, r, fun = 'count')
    plot(mydata.r, col = rev(topo.colors(10)))
    
    # raster with ggplot2
    rasterVis::gplot(mydata.r) +
      geom_tile(aes(fill = value)) +
      scale_fill_gradient(low = 'yellow', high = 'blue') +
      coord_equal()
    

    Retrieve frequencies in the raster

    # Get all data of the raster as a matrix
    m <- as.matrix(mydata.r)
    
    # Retrieve frequencies in more usable way
    # /!\ Be careful, `Var1=1` corresponds to `row=1`
    # you can retrieve real coordinates with `xFromCol` or `yFromRow`
    m2 <- as.data.frame.table(
      as.matrix(mydata.r), 
      base = list(as.character(1:max(xmax(mydata.sp), ymax(mydata.sp))))) %>%
      filter(!is.na(Freq))
    

    There is also a solution with ggplot2::geom_bin2d directly on the original dataframe.

    # ==== Direct with ggplot2 ====
    ggplot(mydata, aes(Position_X, Position_Y)) +
      geom_bin2d(binwidth = c(10, 10)) +
      scale_fill_gradient(low = 'yellow', high = 'blue') +
      coord_equal()
    

    enter image description here