Проблема с наложением кригенных данных на карту Google - PullRequest
Новая
       19

Проблема с наложением кригенных данных на карту Google

0 голосов
/ 22 мая 2018

В приведенном ниже коде я пытаюсь наложить данные kriged meuse на карту Google с помощью ggmap(). Код, кажется, работает нормально вплоть до get_map(),, но застревает в каких-то огромных вычислениях (использование ОЗУ)скачет на 7G на моей машине 8G - будьте осторожны!) при последнем print(ggmap()) выполнении.Что я делаю не так?

ОБНОВЛЕНИЕ: Я думаю, что я выяснил проблему => stat_contour(), похоже, не работает с кривыми данными (все еще не уверен, почему).Когда я заменяю его на функцию binning stat_summary_2d(),, код работает.Я все еще хотел бы показать контурные линии, если это возможно.Может кто-то помочь с этим?

Спасибо 

# transform meuse data to SpatialPointsDataFrame
suppressMessages(library(sp))
data(meuse)
coordinates(meuse) <- ~ x + y
proj4string(meuse) <- CRS("+proj=stere
                          +lat_0=52.15616055555555 +lon_0=5.38763888888889
                          +k=0.999908 +x_0=155000 +y_0=463000
                          +ellps=bessel +units=m +no_defs
                          +towgs84=565.2369,50.0087,465.658,
                          -0.406857330322398,0.350732676542563,-1.8703473836068, 4.0812")

# define sample grid for kriging
set.seed(42)
grid <-spsample(meuse, type = "regular", n = 10000)

# do kriging
suppressMessages(library(automap))
krg <- autoKrige(formula = copper ~ 1,
                 input_data = meuse,
                 new_data = grid)

# extract kriged data
krg_df <- data.frame(krg$krige_output@coords,
                     pred = krg$krige_output@data$var1.pred)
names(krg_df) <- c("x", "y", "pred")

# transform to SpatialPointsDF & assign original (meuse) projection
krg_spdf <- krg_df
coordinates(krg_spdf) <- ~ x + y 
proj4string(krg_spdf) <- proj4string(meuse)

# transform again to longlat coordinates (for overlaying on google map below)
krg_spdf <- spTransform(krg_spdf, CRS("+init=epsg:4326"))
krg_df <- data.frame(krg_spdf@coords, pred = krg_spdf@data$pred)

# get meuse map and overlay kriged data
suppressMessages(library(ggmap))
lon <- range(krg_df$x)
lat <- range(krg_df$y)
meuse_map <- get_map(location = c(lon = mean(lon), lat = mean(lat)),
                     zoom = 13)
print(ggmap(meuse_map, extent = "normal", maprange = F) +
        #### old code (does not work) ####
        #stat_contour(aes(x = x, y = y, z = pred, fill = ..level..),
        #             alpha = 0.5,
        #             color = "gray80",
        #             geom = "polygon",
        #             data = df) +
        #### replaced with stat_summary_2d() ####
        stat_summary_2d(aes(x = x, y = y, z = pred),
                        binwidth = c(0.0005,0.0005),
                        alpha = 0.5,
                        data = krg_df) +
        scale_fill_gradient(low = "yellow",
                            high = "red",
                            name = "Copper") +
        coord_cartesian(xlim = lon, ylim = lat, expand = 0) +
        theme(aspect.ratio = 1))

enter image description here

...