# ЗАКРУЗКА БИБЛИОТЕК
if(!require(dendextend)) install.packages("dendextend")
if(!require(DescTools)) install.packages("DescTools")
if(!require(dismo)) install.packages("dismo")
if(!require(dplyr)) install.packages("dplyr")
if(!require(geosphere)) install.packages("geosphere")
if(!require(ggalt)) install.packages("ggalt")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(ggthemes)) install.packages("ggthemes")
if(!require(GISTools)) install.packages("GISTools")
if(!require(Hmisc)) install.packages("Hmisc")
if(!require(raster)) install.packages("raster")
if(!require(rgdal)) install.packages("rgdal")
if(!require(rgeos)) install.packages("rgeos")
if(!require(sp)) install.packages("sp")
if(!require(spatstat)) install.packages("spatstat")
if(!require(spgwr)) install.packages("spgwr")
if(!require(tibble)) install.packages("tibble")
if(!require(utils)) install.packages("utils")

library(dendextend)
library(DescTools)
library(dismo)
library(dplyr)
library(geosphere)
library(ggalt)
library(ggplot2)
library(ggthemes)
library(GISTools)
library(Hmisc)
library(raster)
library(rgdal)
library(rgeos)
library(sp)
library(spatstat)
library(spgwr)
library(stats)
library(tibble)
library(utils)

# 1. ЗАГРУЗКА ДАННЫХ
# a) загрузка карты и данных
#datMOSPb - данные муниципальных округов
datMOSPb <- csv.get("http://lornii.ru/resources/lib/R/MOSPb.csv", 
                    vnames=1, # первая строка в csv - переменные
                    #labels=2, # вторая метки
                    skip=1, # две верхние строки пропускаем - там названия переменных и метки
                    sep=";") # разделитель в Excel для csv - ";"
datMOSPb<-as.data.frame(datMOSPb)
#stpat - данные о пролеченных больных в стационаре
stpat <- csv.get("http://lornii.ru/resources/lib/R/spatxy2018.csv", 
                 vnames=1, # первая строка в csv - переменные
                 labels=2, # вторая метки
                 skip=2, # две верхние строки пропускаем - там названия переменных и метки
                 sep=";",
                 stringsAsFactors=FALSE) # разделитель в Excel для csv - ";"
stpat<-as.data.frame(stpat)
stpatsp<-as.data.frame(stpat)
# geoMOSPb - полигоны муниципальных образований (МО)
geoMOSPb <- readOGR("http://lornii.ru/resources/lib/R/datMOSPb.geojson", "datMOSPb", require_geomType="wkbPolygon", use_iconv=TRUE, encoding="UTF-8", stringsAsFactors = FALSE)
# преобразование  полигонов МО для отображения на диаграммах ggplot2
geoMOSPb_fort <- fortify(geoMOSPb, region = "id") # to dataframe обязательно указать region

# b) Субсеттинг петербуржцев с сенсоневральной тугоухостью (СНТ)
spSPb<-subset(stpat, AREA=="Санкт-Петербург")
patternsTU <- c("H90.3", "H90.4", "H90.5", "H91.1", "H91.2", "H91.8") #сенсоневральная тугоухость (СНТ)
spSPbTU <- filter(spSPb, grepl(paste(patternsTU, collapse="|"), MKB))

# СОЗДАНИЕ КАРТОГРАММ 
#для Рис. 1: а) численности населений муниципальных образований СПб, б) частоты госпитализации (отн. значение, 0/00) в НИИ ЛОР на 1 тыс. жителей МО СПб и в) общего числа госпитализированных (абс. значение) в НИИ ЛОР из МО СПб
#ПОСЧЕТ КОЛИЧЕСТВА ПАЦИЕНТОВ В ПОЛИГОНАХ
# Создаем набор координат загрубленных широты и долготы места жительства пациентов для подсчета их количества в полигонах и преобразуем их в объект "SpatialPointsDataFrame"
myvars<-c("lon","lat")
xyTU <- spSPbTU[myvars]
xyTU$lat<-as.numeric(xyTU$lat)
xyTU$lon<-as.numeric(xyTU$lon)
spspSPbTU <- SpatialPointsDataFrame(coords = xyTU, data = spSPbTU,
                                    proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

# Подсчитываем  число пациентов в полигонах МО, создаем переменную с их количеством 'countpatTU' и вычисляем распространенность случаев госпитализации  больных с СНТ из МО 
# Подсчитываем число пациентов СНТ в полигонах МО
countTUmo<-poly.counts(pts = spspSPbTU, polys = geoMOSPb) 
npatTUinmo<-setNames(countTUmo, geoMOSPb@data$id)
countpatTUmo <- stack(npatTUinmo) 
countpatTUmo$ind<-as.character(countpatTUmo$ind)
datMOSPb$relOSM<-as.character(datMOSPb$relOSM)
sppop<-left_join(datMOSPb, countpatTUmo, by = c("relOSM" = "ind"))
names(sppop)[names(sppop) == 'values'] <- 'countpatTU'
sppop$countpatTU<-as.numeric(sppop$countpatTU)
sppop$prevalenceTU<-(sppop$countpatTU/sppop$population)*1000
sppop$prevalenceTU<-as.numeric(sppop$prevalenceTU)

# Находим центроиды муниципальных округов, к которым будем привязывать данные распространенности  болезней и т.п.
centroids <- as.data.frame(gCentroid(geoMOSPb, byid = TRUE, id = geoMOSPb@data$id))
colnames(centroids) <- c("lon", "lat") 
# Объединяем координаты центроидов муниципальных округов с основной таблицей МО
sppop<-left_join(sppop, tibble::rownames_to_column(centroids), by = c("relOSM" = "rowname"))
geoMOSPb<-merge(geoMOSPb, sppop[,c("relOSM", "population","countpatTU", "prevalenceTU", "lon","lat")], by.x="id", by.y= "relOSM",all.x=TRUE)
geoMOSPb@data<-geoMOSPb@data[,c("id","name","place", "population.y","countpatTU",  "prevalenceTU",  "lon","lat")]
names(geoMOSPb@data)[names(geoMOSPb@data) == 'population.y'] <- 'population' # поменял имя столбца№

# Создаем картограммы Рис. 1
# Рис. 1 а)
ggplot() +
  geom_cartogram(data = geoMOSPb_fort, aes(x = long, y = lat, map_id = id), 
                 map = geoMOSPb_fort) +
  geom_cartogram(data = sppop, aes(fill = population/1000, map_id = relOSM),
                 map = geoMOSPb_fort, color = "black", size = 0.3) +
  scale_fill_gradientn(name="Численность \nжителей (тыс. чел.).", colours = rev(brewer.pal(7, "Spectral"))) +
  coord_map() +
  theme_map()
# Рис. 1 б)
ggplot() +
  geom_cartogram(data = geoMOSPb_fort, aes(x = long, y = lat, map_id = id), 
                 map = geoMOSPb_fort) +
  geom_cartogram(data = sppop, aes(fill = prevalenceTU, map_id = relOSM),
                 map = geoMOSPb_fort, color = "black", size = 0.3) +
  scale_fill_gradientn(name="Число госпитализированных \nc СНТ на 1000 жителей", colours = rev(brewer.pal(7, "Spectral"))) +
  coord_map() +
  theme_map()
# Рис. 1 в)
ggplot() +
  geom_cartogram(data = geoMOSPb_fort, aes(x = long, y = lat, map_id = id), 
                 map = geoMOSPb_fort) +
  geom_cartogram(data = sppop, aes(fill = countpatTU, map_id = relOSM),
                 map = geoMOSPb_fort, color = "black", size = 0.3) +
  scale_fill_gradientn(name="Число госпитализированных", colours = rev(brewer.pal(7, "Spectral"))) +
  coord_map() +
  theme_map()

#3. ТОЧЕЧНАЯ ДИАГРАММА ПАЦИЕНТОВ С СНТ И ОЦЕНКА ПЛОТНОСТИ ВЕРОЯТНОСТИ ИХ ПОЯВЛЕНИЯ (Рис. 2)
spSPbTU$lat<-as.numeric(spSPbTU$lat) # из "labelled"  "character"
spSPbTU$lon<-as.numeric(spSPbTU$lon)
g1<-ggplot() +
  geom_polygon(data = geoMOSPb_fort, aes( x = long, y = lat, group = group), fill="white", color="grey") +
  theme_void() +
  coord_map()
g1 + 
  geom_density_2d(aes(x = lon, y = lat),data = spSPbTU)+
  geom_point(aes(x = lon, y = lat),size = 2, colour="red", shape=24, show.legend = T,data = spSPbTU)+
  theme(legend.position = 'none')

#4. ВЫЧИСЛЕНИЕ K и L ФУНКЦИЙ Рис. 3 
spSPbTU.ppp <- as(SpatialPoints(spspSPbTU), "ppp")
Kloh<-lohboot(spSPbTU.ppp, Kest) # К-функция
plot(Kloh)
Lloh<-lohboot(spSPbTU.ppp, Lest) # б) L-функция
plot(Lloh)

#5. РАСЧЕТ МОДЕЛИ ГЕОГРАФИЧЕСКОЙ ВЗВЕШЕННОЙ РЕГРЕССИИ Рис. 4 
GWRbandwidth <- gwr.sel(countpatTU ~ population, data=geoMOSPb@data, coords=cbind(geoMOSPb@data$lon,geoMOSPb@data$lat),adapt=T) 
gwr.model = gwr(countpatTU ~ population, data=geoMOSPb@data, coords=cbind(geoMOSPb@data$lon,geoMOSPb@data$lat),adapt=GWRbandwidth, hatmatrix=TRUE, se.fit=TRUE) 
gwr.model

#6 КАРТОГРАММА КОЭФФИЦИЕНТОВ GWR С ПРИВЯЗКОЙ К МО СПб Рис. 5.
results<-as.data.frame(gwr.model$SDF)
head(results)
geoMOSPb@data$coefpopulation<-results$population
ggplot() +
  geom_cartogram(data = geoMOSPb_fort, aes(x = long, y = lat, map_id = id), 
                 map = geoMOSPb_fort) +
  geom_cartogram(data = geoMOSPb@data, aes(fill = coefpopulation, map_id = id),
                 map = geoMOSPb_fort, color = "black", size = 0.3) +
  scale_fill_gradientn(name="Коэффициент \n регрессии GWR", colours = rev(brewer.pal(7, "Spectral"))) +
  coord_map() +
  theme_map()

#7 РЕЗУЛЬТАТЫ ИЕРАРХИЧЕСКОГО КЛАСТЕРНОГО АНАЛИЗА Рис. 6

mdist <- distm(spspSPbTU)
hc <- hclust(as.dist(mdist), method="complete")
d=500
k=4
spspSPbTU$clust <- cutree(hc, h=d, k = k)
spspSPbTU@bbox[] <- as.matrix(extend(extent(spspSPbTU),0.001))
# get the centroid coords for each cluster
cent <- matrix(ncol=2, nrow=max(spspSPbTU$clust))
for (i in 1:max(spspSPbTU$clust))
  cent[i,] <- gCentroid(subset(spspSPbTU, clust == i))@coords # gCentroid from the rgeos package
ci <- circles(cent, d=d, lonlat=T) #library(dismo)
clcol<-rainbow(k)[factor(spspSPbTU$clust)]
plot(geoMOSPb,border="grey")
plot(spspSPbTU, col=clcol,  add=T)
plot(ci@polygons, axes=T, add=T, col="black")
#legend("bottomleft", legend = c("1","2","3","4"),  
#       pch = "+", col = c("red","blue", "green","yellow"),
#       bty = "n", title = "Кластеры")