Dwa dni temu, w tym wpisie przedstawialiśmy grafiki przygotowane przez Michała Kurtysa w ramach wakacyjnego projektu.
Pan Michał przygotował też krótki tutorial wyjaśniający jak samodzielnie zrobić takie wykresy w R. Poniżej wklejamy ten tutorial. Jest on moim zdaniem bardzo ciekawy i porusza wiele ciekawych technicznych problemów z analizą danych przestrzennych.
Dane
Informacje o kolejkach zapisane są w plikach Excelowych o rozszerzeniu xls.
Na każde województwo przypada takich plików kilka. Nas teraz będzie interesował tylko plik dotyczący świadczeń specjalistycznych, zawierający w nazwie skrót AOS.
W sumie jest więc 16 plików (po jednym na województwo), które wyglądają mniej-więcej w ten sposób.
01_AOS_31072013.xls
02_AOS_31072013.xls
...
16_AOS_31072013.xls
Województwa są ułożone alfabetycznie – najmniejszy numer odpowiada województwu dolnośląskiemu, a największy zachodniopomorskiemu. Aby powiązać numer z nazwą województwa przygotowałem plik województwa.csv, w którym są one wypisane w należytej kolejności.
wojwództwa.csv
WOJ. DOLNOŚLĄSKIE
WOJ. KUJAWSKO-POMORSKIE
WOJ. LUBELSKIE
WOJ. LUBUSKIE
WOJ. ŁÓDZKIE
WOJ. MAŁOPOLSKIE
WOJ. MAZOWIECKIE
WOJ. OPOLSKIE
WOJ. PODKARPACKIE
WOJ. PODLASKIE
WOJ. POMORSKIE
WOJ. ŚLĄSKIE
WOJ. ŚWIĘTOKRZYSKIE
WOJ. WARMIŃSKO-MAZURSKIE
WOJ. WIELKOPOLSKIE
WOJ. ZACHODNIOPOMORSKIE
Dane o granicach administracyjnych Polski pobrałem z Geoportalu. Ze względów licencyjnych początkowo chciałem wykorzystać OpenStreetMap. Gotowe pliki shapefile oferuje m.in. geofabrik/cloudmade. Niestety, jeżeli chodzi o Polskę nie miałem tam czego szukać – wszystkie mapy były zdeformowane.
Mapy z Geoportalu nie są dostępne w formacie shapefile. Na szczęście dzięki instrukcjom we wpisie Pawła Wiechuckiego przedstawianego na tym blogu ich własnoręczne stworzenie nie było trudne.
Biblioteki
Potrzebujemy następujących bibliotek:
** maptools – funkcja readShapePoly pozwala na wczytanie pliku ShapeFile.
** sp – do manipulacji i wyświetlania danych kartograficznych.
** rgeos, rgdal – funkcja spTransform, pozwala zmienić układ współrzędnych
** FNN – funkcja get.knnx – do regresji
** SmarterPoland – funkcję getGoogleMapsAddress ułatwi pobranie współrzędnych punktu o określonym adresie
** Cairo – do produkcji wykresów
** XLConnect – obsługa plików excela
Kod R
Zacznijmy od wczytania bibliotek.
1
2
3
4
5
6
7
8
| library(maptools)
library(sp)
library(rgeos)
library(SmarterPoland)
library(FNN)
library(rgdal)
library(Cairo)
library(XLConnect) |
library(maptools)
library(sp)
library(rgeos)
library(SmarterPoland)
library(FNN)
library(rgdal)
library(Cairo)
library(XLConnect)
Następnie zmiana katalogu roboczego. Wczytanie listy województw.
Z każdego pliku „aos” odczytujemy dane, które zaczynają się w 3 wierszu(razem z nagłówkiem) i mają 9 kolumn.
Postawnowiłem zmienić nazwy kolumn, gdyż oryginalne były niesłychanie długie.
Dodatkowo tworzymy dodatkowe kolumny – ID i nazwę województwa.
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
| 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") |
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")
Pacjenci są podzieleni przez NFZ na dwie grupy – „przypadek stabilny” i „przypadek ostry”.
Graficznie przedstawiać będziemy wyłącznie dane o przypadkach stabilnych.
Dalej dzielimy adres na części składowe – nazwa miejscowości i ulica razem z numerem numer lokalu.
36
37
38
39
| 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) |
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)
Pobieramy współrzędne wyłącznie miejscowości w której znajduje się poradnia.
Jest to znacznie szybsze, a z pewnością nie potrzebujemy większej dokładności.
Po za tym Google ogranicza darmowy dostęp do usługi do 2500 zapytań dziennie.
Konstruujemy więc tabelę zawierającą nazwę miejscowości w jednej kolumnie, a w drugiej województwo w którym się ona znajduje.
Dodanie województwa pomaga uściślić zapytanie – istnieją przecież miejscowości, które noszą tę samą nazwę.
Funkcja unique eliminuje zduplikowane wiersze.
W pętli tworzymy nową zmienną region, w której zapisujemy jedynie nazwę województwa bez skrótu „Woj.” na początku.
Google Geocoding raczej nie trawi tego przedrostka.
Województwo łączymy z nazwą miejscowości i przekazujemy do argumentu city funkcji getGoogleMapAddress.
Może to wyglądać to dziwnie, ale zapytanie i tak zostanie skonstruowane poprawnie.
40
41
42
43
44
45
46
47
48
49
50
51
52
| 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_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]
}
Zawczasu warto przygotować sobie współrzędne w układzie, który jest bardziej przyjazny do zadania:
Tworzymy zmienną typu SpatialPoints i określamy ich układ współrzędnych.
Współrzędne, które pobraliśmy są w układzie odniesienia WGS 84, znanym także jako EPSG:4326
Zmieniamy układ współrzędnych przy pomocy funkcji spTransform.
Ostatecznie łączymy dwie tabele przy pomocy funkcji merge.
Wreszcie mamy wszystkie interesujące nas wartości.
53
54
55
56
57
58
59
60
| 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) |
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)
Wczytanie pliku shapefile, który zawiera granice administracyjne województw.
Jest już w prawidłowym układzie współrzędnych – wystarczy tylko go określić.
61
62
| wojewodztwa.shp = readShapePoly("geoportal\\woj.shp")
proj4string(wojewodztwa.shp)=CRS("+init=epsg:2180") |
wojewodztwa.shp = readShapePoly("geoportal\\woj.shp")
proj4string(wojewodztwa.shp)=CRS("+init=epsg:2180")
Funkcja generate_grid posłuży do wygenerowania kraty, która będzie pokrywać powierzchnię określonego pola.
Na obrazku widać punkty kraty. W istocie są to punkty, które będą służyły jako punkty startowe linii.
Pierwszy argument funkcji to ilość komórek kraty w osi X.
Komórka kraty ma być kwadratem, więc ich ilość w osi Y nie będzie przekazywana do funkcji.
Następnie wyliczamy topologię kraty. Przekazujemy:
** współrzędne dolnego, lewego rogu
** wymiary komórek
** ilość komórek
I na koniec zwracamy obiekt klasy SpatialGrid.
63
64
65
66
67
68
69
70
71
72
73
| 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)
} |
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)
}
Jak widać na obrazku, część punktów kraty znajduje się poza granicami kraju.
Funkcja which_cells_inside zwraca indeksy komórek kraty, które znajdują się wewnątrz obiektu SpatialPolygons.
Wykorzystamy do tego funkcję over. Dla argumentów SpatialPoints i SpatialPolygons funkcja zwraca wektor o długości równej ilości punktów.
Wartości wektora mówią o tym, wewnątrz którego wielokąta znajduje się punkt.
Tych mamy 16 – odpowiadają województwom.
Jeżeli jest poza jakimkolwiek wielokątem składowym to danemu punktowi odpowiadać będzie wartość NA.
Szczegółowo wielokąt składowy to obiekt klasy Polygons:
74
75
76
77
78
79
80
81
82
83
84
85
| > 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 |
> 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
86
87
88
89
90
91
92
93
| 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)
} |
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)
}
Generowanie linii. W argumentach znajduje się zmienna cell_distance, która mówi co ile komórek kraty będzie znajdować się początek linii.
Gdy pisałem kod uznałem, że być może takie rozwiązanie się przyda. Teraz jestem pewien, że bardziej elegancko byłoby generować mniejszą kratę.
Co z resztą każdy zauważy.
Zwracamy obiekt typu SpatialLines. Strukura obiektów (Spatial)Line(s) jest podobna do Polygons.
Przypominam, że linie rysujemy od punktu kraty do najbliższej przychodni.
Musimy więc wiedzieć, która jest najbliższa.
Przekazujemy więc macierz knn_indices w której będzie to zapisane.
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
| # 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)
} |
# 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)
}
Odpalamy funkcje.
W pętli przejdziemy po kolei po ka¿dym typie poradni.
Wewnątrz niej wybierzemy interesujące nas wiersze w tabeli poradnie.
Dzięki bibliotece Cairo łatwo zapiszemy wykres w wysokiej rozdzielczości.
get.knnx jest funkcją obsługującą regresję k-sąsiadów.
W pierwszym argumencie jest zbiór danych, w drugim zbiór zapytań.
W tej sytuacji interesuje nas wyłącznie najbliższy sąsiad, więc k wynosi 1.
Zwraca dwie macierze – w pierwszej zapisuje indeksy najbliższych sąsiadów, a w drugiej odległości.
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
| #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()
} |
#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()
}
Hidden track. Tworzenie map „kolorowych”. Tworzymy gęstszą kratę.
Następnie uzyskujemy punkty kraty i zmieniamy ich układ współrzędnych do WGS 84.
Jest to niezbędne – dla układu współrzędnych epsg:2180 funkcja spDistsN1 nie zwracała odległości w kilometrach.
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
| 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()
} |
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()
}
Było to mój pierwszy ,,większy” projekt w R. Z pewnością wiele rozwiązań nie jest najlepszych, ale mimo wszystko mam nadzieję, że ten wpis będzie pomocny.
