require(maps) #if required install with install.packages("maps") require(mapdata) #If required install with install.packages("mapdata") distxy <- function(coorda, coordb) #This is our custom function to calculate earth distances { dist <- acos(sin(pi/180*coorda[1])*sin(pi/180*coordb[1]) + cos(pi/180*coorda[1])* cos(pi/180*coordb[1])*cos(pi/180*coordb[2]-pi/180*coorda[2]))*6370 return(dist) } #quartz() # Use windows() when running R on a windows machine #pdf("Graph.pdf") data <- read.csv("koord.csv") # Read required data lim <- c(-126,-66) # Set plotting border lim2 <- c(23, 55) # Set iteration start point xstart <- 45 ystart <- -110 number_it <- 20 iteration_colors <- rainbow(number_it, alpha = 1.0, start=0.33, end=0.7) data2 <- data.frame(x=data[,3], y=data[,4]) xstart_o <- 45 ystart_o <- -110 # plot map and customers it <- data.frame(x=0, y=0) costs <- c(NULL) # optimization iterations for (i in 1:number_it) { distance <- apply(data2, 1, distxy, coordb=c(xstart,ystart)) # calculate current distances nenner <- data[,5]/distance # denominator zahler_x <- data[,5]*data2[,1] / distance #nominator_y zahler_y <- data[,5]*data2[,2] / distance #nominator_y xstart <- sum(zahler_x) / sum(nenner) ystart <- sum(zahler_y) / sum(nenner) it[i,1] <- xstart #save current x_coordinate it[i,2] <- ystart #save current y_coordinate costs<-append(costs, sum(data[,5]*distance)) #calculate costs of current solution } # Now plot the results map("state", lwd=0.3, mar=c(0,0,0,0), col="#ffe895", border=1, interior=TRUE,fill=TRUE, bg="#93c3f2", resolution=0, , xlim = lim, ylim= lim2,) points(data[,4], data[,3], col="darkred", pch=21,bg="red", cex = sqrt(data[,5]/3)*0.1) points(ystart_o, xstart_o, col="darkblue", pch=16, cex=1) # Startpunkt zeichnen points(it[,2], it[,1], col = iteration_colors[number_it:1], pch=16) #dev.off()