library(maptools)
library(sp)
library(rgeos)
library(SmarterPoland)
library(FNN)
library(rgdal)
library(Cairo)
library(XLConnect)

setwd("C:\\Users\\mic\\Documents\\eurostat\\nfzqueue\\new")
wojewodztwa = read.csv("wojewodztwa.csv", header=FALSE)
wojewodztwa$V1 = as.character(wojewodztwa$V1)
 
poradnie = data.frame()
for( i in 1:16) {
	filename = sprintf("%02d_AOS_31072013.xls", i)
	print(filename)
	wb = loadWorkbook( filename, create = FALSE)
	dane = readWorksheet(wb,sheet="Zestawienie", startRow = 3, startCol=0, endCol=8)
	dane$IDWojewodztwa = i
	dane$Wojewodztwo = wojewodztwa$V1[i]
	poradnie = rbind(poradnie, dane)
}
 
#zmiana nazwy kolumn
names(poradnie) = c(
"Nazwa.komorki.realizujacej",
"Kategoria",
"Nazwa",
"Nazwa.komorki",
"Adres",
"Liczba.oczekujacych",
"Liczba.skreslonych",
"Sredni.czas",
"ID.wojewodztwa",
"Wojewodztwo")

poradnie = poradnie[which(poradnie$Kategoria == "przypadek stabilny"), ]
poradnie$Miejscowosc = sapply(strsplit(poradnie$Adres, "\n"), function(x) x[1] )
poradnie$Ulica = sapply(strsplit(poradnie$Adres, "\n"), function(x) x[2] )
poradnie$Sredni.czas = as.numeric(poradnie$Sredni.czas)

places_unique = data.frame( poradnie$Miejscowosc, poradnie$Wojewodztwo, check.names=FALSE)
places_unique = unique( places_unique )
names(places_unique) = c("Miejscowosc", "Wojewodztwo")
places_unique$Wojewodztwo = as.character(places_unique$Wojewodztwo)
places_unique$Miejscowosc = as.character(places_unique$Miejscowosc)
 
for( i in 1:nrow(places_unique) ) {
	region = strsplit( places_unique$Wojewodztwo[i], " " )[[1]][2]
	temp_coords = getGoogleMapsAddress(street="",city=paste(places_unique$Miejscowosc[i],",",region))
 
	places_unique$lat_wgs84[i] = temp_coords[1]
	places_unique$lng_wgs84[i] = temp_coords[2]
}

places_sp = SpatialPoints( cbind(places_unique$lng_wgs84, places_unique$lat_wgs84) )
proj4string(places_sp) = CRS("+init=epsg:4326")
places_sp = spTransform( places_sp, CRS("+init=epsg:2180"))
 
places_unique$lng = coordinates(places_sp)[,1]
places_unique$lat = coordinates(places_sp)[,2]
 
poradnie = merge(poradnie, places_unique)

wojewodztwa.shp = readShapePoly("geoportal\\woj.shp")
proj4string(wojewodztwa.shp)=CRS("+init=epsg:2180")

generate_grid = function( cell_n_x=1000, sp_polygons) {
	bbox_min = bbox(sp_polygons)[,1]
	bbox_max = bbox(sp_polygons)[,2]
	cell_width=((bbox_max-bbox_min)/cell_n_x)[1]
	cell_height=cell_width
	cell_n_y = ceiling(((bbox_max-bbox_min)/cell_height)[2])
 
	gt = GridTopology( bbox_min, c(cell_width,cell_height), c(cell_n_x,cell_n_y)  )
	sg = SpatialGrid(gt, CRS("+init=epsg:2180"))
	return(sg)
}

> class(wojewodztwa.shp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> class(wojewodztwa.shp@polygons)
[1] "list"
> class(wojewodztwa.shp@polygons[[1]])
[1] "Polygons"
attr(,"package")
[1] "sp"
> length(wojewodztwa.shp@polygons)
[1] 16

which_cells_inside = function( sp_grid, sp_polygons) {
	grid_coords = coordinates(sp_grid)
	grid_points = SpatialPoints(grid_coords, CRS( proj4string(sp_grid)) )
	inside = over( grid_points, sp_polygons)
	proj4string(grid_points) = proj4string(wojewodztwa.shp)
	inside = which( inside$ID1 > 0 )
	return(inside)
}

# sp_points to obiekt typu SpatialPoints, zawierający listę przychodni.
# cell_subset ogranicza kratę, tak by linie były prowadzone tylko z wewnątrz kraju. 
 
generate_lines = function( knn_indices, sp_grid, sp_points, cell_distance=1, cell_subset=NULL ) {
	cell_n_x = sg@grid@cells.dim[1]
	cell_n_y = sg@grid@cells.dim[2]
	#macierz ktora zawiera indeksy komorek w ktorych zaczynaja sie linie
	line_anchors = matrix(nrow=floor(cell_n_y/cell_distance), ncol=floor(cell_n_x/cell_distance))
	for( i in 1:nrow(line_anchors)) {
		line_anchors[i,] = seq(cell_distance, cell_n_x, by= cell_distance) +
						  rep(cell_n_x*cell_distance*i,floor(cell_n_x/cell_distance))
	}
	#tylko te linie ktorych poczatek jest w wybranych komorkach
	#na przyklad wewnatrz kraju
	line_anchors = intersect(line_anchors, cell_subset)
 
	#konstrukcja obiektow biblioteki sp
	line_starts = coordinates(sp_grid)[line_anchors,]
	closest_medic_idx = knn_indices[line_anchors,1]
 
	line_stops = coordinates(sp_points)[closest_medic_idx,]
	line_list = sapply( 1:nrow(line_starts), function(x) Line( rbind( line_starts[x,], line_stops[x,])) )
	line_obj = Lines(line_list, ID="a")
	sp_lines_obj = SpatialLines( list(line_obj) )
 
	return(sp_lines_obj)
}

#wczytanie danych, etc..
#kod
sg = generate_grid( cell_n_x=100, wojewodztwa.shp)
inside_poland = which_cells_inside( sg, wojewodztwa.shp)
 
for( specialist_type in unique( (poradnie$Nazwa.komorki.realizujacej) ) ) {
	print(specialist_type)
	flush.console()
	specialist = poradnie[which(poradnie$Nazwa.komorki.realizujacej == specialist_type),]
	specialist_sp = SpatialPoints( cbind( specialist$lng, specialist$lat) ,CRS("+init=epsg:2180") )
 
	knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = 1 )
 
	sp_lines_obj = generate_lines( knn_res[[1]], sg, specialist_sp, cell_distance=1, cell_subset=inside_poland )
 
	Cairo(900, 700, file=paste(specialist_type,".png",sep=""), type="png", bg="white")
	plot(wojewodztwa.shp)
	plot(sp_lines_obj,add=T)
	#plot(specialist_sp,col="blue",lwd=2, add=T)
	title(main=specialist_type)
	dev.off()
}

sg = generate_grid( cell_n_x=500, wojewodztwa.shp)
inside_poland = which_cells_inside( sg, wojewodztwa.shp)
 
grid_points_wgs84 = SpatialPoints( coordinates(sg), CRS(proj4string(sg)))
grid_points_wgs84 = spTransform(grid_points_wgs84, CRS("+init=epsg:4326"))
gp_coords = coordinates(grid_points_wgs84)
 
for( specialist_type in unique( (poradnie$Nazwa.komorki.realizujacej) ) ) {
	print(specialist_type)
	flush.console()
	specialist = poradnie[which(poradnie$Nazwa.komorki.realizujacej == specialist_type),]
	specialist_sp = SpatialPoints( cbind( specialist$lng, specialist$lat) ,CRS("+init=epsg:2180") )
	specialist_coords_wgs84 = cbind( specialist$lng_wgs84, specialist$lat_wgs84)
 
	#knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = ifelse(nrow(specialist)>10, 10, nrow(specialist) ))
	knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = 1)
 
	knn_indices = knn_res[[1]]
	km_distances =  matrix( nrow=nrow( knn_indices ), ncol=ncol(knn_indices))
	closest_coords = apply( knn_indices, c(1,2), function(x) specialist_coords_wgs84[x,])
	closest_coords = matrix( closest_coords, ncol = 2, byrow=TRUE ) #2*ncol(knn_indices)
 
	for (rowID in 1:nrow( knn_indices ) ) {
		pts = matrix(closest_coords[rowID,],ncol=2)
 
		p = gp_coords[rowID,]
		km_distances[rowID,] = spDistsN1(pts, p, TRUE)
	}
	closest_distance = km_distances[,1]
	sg_df = SpatialGridDataFrame( sg, data.frame(closest_distance))
 
	#obszar położony poza granicami kraju ma kolor biały
	sg_df@data$closest_distance[-inside_poland]=0
	#80 jest górną granicą wartości na wykresie
	sg_df@data$closest_distance[ sg_df@data$closest_distance>80] = 80
 
	Cairo(900, 700, file=paste(specialist_type,"_c",".png",sep=""), type="png", bg="white")
	print(spplot( sg_df, panel = function(x,y,...){
	panel.gridplot(x,y,...,  )
	sp.polygons( wojewodztwa.shp )
	sp.points( specialist_sp, pch=".", col="black" )
	}, at=c(1:80), col.regions = rev(heat.colors(100)), main=specialist_type))
	dev.off()
 
}