Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Reproducible example for my SO question on how to fit 2 sf seamlessly -- https://stackoverflow.com/questions/48360446
# install dev version of ggplot2
devtools::dev_mode()
devtools::install_github("tidyverse/ggplot2")
library(tidyverse)
library(sf)
library(rmapshaper)
library(ggthemes)
# load data
source(file = url("https://gist.githubusercontent.com/ikashnitsky/4b92f6b9f4bcbd8b2190fb0796fd1ec0/raw/1e281b7bb8ec74c9c9989fe50a87b6021ddbad03/minimal-data.R"))
# test how good they fit together
ggplot() +
geom_sf(data = REG, color = "black", size = .2, fill = NA) +
geom_sf(data = NEI, color = "red", size = .2, fill = NA)+
coord_sf(datum = NA)+
theme_map()
ggsave("test-1.pdf", width = 12, height = 10)
# simplify
REGs <- REG %>% ms_simplify(keep = .5, keep_shapes = TRUE)
NEIs <- NEI %>% ms_simplify(keep = .5, keep_shapes = TRUE)
ggplot() +
geom_sf(data = REGs, color = "black", size = .2, fill = NA) +
geom_sf(data = NEIs, color = "red", size = .2, fill = NA)+
coord_sf(datum = NA)+
theme_map()
ggsave("test-2.pdf", width = 12, height = 10)
################################################################################
#
# ikashnitsky.github.io 2018-01-20
# Reproducible example to the SO question
# https://stackoverflow.com/questions/48360446
# How to fit 2 sf seamlessly?
# Ilya Kashnitsky, ilya.kashnitsky@gmail.com
#
################################################################################
# install dev version of ggplot2
devtools::dev_mode()
devtools::install_github("tidyverse/ggplot2")
library(tidyverse) # R 2.0
library(lubridate) # deal with times
library(janitor) # auto clean dataframes
library(sf) # GIS in R 2.0
library(rmapshaper) # to simplify polygons nicely
library(rgdal) # to read the shapefile not as sf
library(ggthemes) # theme_map()
# geodata will be stored in a directory "geodata"
ifelse(!dir.exists('geodata'),
dir.create('geodata'),
paste("Directory already exists"))
# Download geodata
# Eurostat official shapefiles for regions
# http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units
# download zipped shapefile (1.9 MB)
f <- tempfile()
download.file("http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2013_20M_SH.zip", destfile = f)
unzip(f, exdir = "geodata/.")
# remote areas to remove (NUTS-2)
remote <- c(paste0('ES',c(63,64,70)),paste('FRA',1:5,sep=''),'PT20','PT30')
# read and clean the shapefile
gd_nuts <- st_read("geodata/NUTS_2013_20M_SH/data/.", "NUTS_RG_20M_2013") %>%
clean_names() %>%
# re-project using EPSG:3035, an appropreate projection for Europe
st_transform(3035) %>%
# filter out remote areas
filter(!str_sub(nuts_id, 1, 4) %in% remote,
!str_sub(nuts_id, 1, 2) == "AL") %>%
mutate(nuts_id = nuts_id %>% as.character())
# filter NUTS-3
gd_nuts3 <- gd_nuts %>% filter(stat_levl==3)
# aggregate NUTS-2 from NUTS-3
gd_nuts2 <- gd_nuts3 %>%
mutate(nuts_id = nuts_id %>% str_sub(1, 4)) %>%
ms_dissolve(field = "nuts_id")
# aggregate NUTS-0 from NUTS-2
gd_nuts0 <- gd_nuts2 %>%
mutate(nuts_id = nuts_id %>% str_sub(1, 2)) %>%
ms_dissolve(field = "nuts_id")
# test plot NUTS-2 with a random color
gd_nuts2 %>%
mutate(rcol = rgb(runif(310), runif(310), runif(310))) %>%
ggplot()+
geom_sf(aes(fill = rcol), color = NA)+
scale_fill_identity() +
theme_map()
# borders -----------------------------------------------------------------
gd_borders <- gd_nuts0 %>%
ms_innerlines() %>%
as_tibble() %>%
st_as_sf()
ggplot() +
geom_sf(data = gd_nuts0, color = "brown") +
geom_sf(data = gd_borders, color = "yellow", size = 1)
# get neighbouring countries ----------------------------------------------
# let's add neighbouring countries
f <- tempfile()
download.file("http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/CNTR_2010_20M_SH.zip", destfile = f)
unzip(f, exdir = "geodata/.")
# filter without European countries - to reduce the final size of geodata
neigh_subset <- c("IM", "FO", "GI", "AD", "MC", "VA", "SM", "BA", "AL", "RS",
"MD", "UA", "BY", "RU", "EG", "LY", "TN", "DZ", "MA", "GG",
"JE", "KZ", "AM", "GE", "AZ", "SY", "IQ", "IR", "IL", "JO",
"PS", "SA", "LB", "MN", "LY", "EG")
# cut off everything behind the map borders
rect <- rgeos::readWKT(
"POLYGON((25e5 13.5e5,
75e5 13.5e5,
75e5 55.5e5,
25e5 55.5e5,
25e5 13.5e5))"
) %>%
st_as_sf() %>%
`st_crs<-`(3035)
# read in the world shapefile as sf object and do all the cleaning
gd_neigh <- st_read("geodata/CNTR_2010_20M_SH/CNTR_2010_20M_SH/Data/.",
"CNTR_RG_20M_2010") %>%
clean_names() %>%
filter(cntr_id %in% neigh_subset) %>%
st_transform(3035) %>%
ms_dissolve() %>%
st_intersection(rect) # cut behind the map
# test how god they fit together
ggplot() +
geom_sf(data = gd_nuts2, color = "black", size = 1, fill = NA) +
geom_sf(data = gd_neigh, color = "red", size = .2, fill = NA)+
theme_map()
ggsave("test-1.pdf", width = 25, height = 20)
# simplify polygons -------------------------------------------------------
# merge neighbours trick --------------------------------------------------
# treat the dissolved polygon for neighbouring countries as a single NUTS-3
gd_trick <- gd_neigh %>%
transmute(nuts_id = "NEIGH",
stat_levl = 3L,
geometry) %>%
rbind(gd_nuts3 %>% select(nuts_id, stat_levl, geometry))
# simp NUTS-3
gd_trick_s <- gd_trick %>% ms_simplify(keep = .5, keep_shapes = TRUE)
gd_n3s_nei <- gd_trick_s %>% filter(nuts_id == "NEIGH")
gd_n3s <- gd_trick_s %>% filter(!nuts_id == "NEIGH")
ggplot() +
geom_sf(data = gd_n3s, size = .1)+
geom_sf(data = gd_n3s_nei, size = .1, color = "red")+
theme_map()
ggsave("test-2.pdf", width = 25, height = 20)
REG = structure(list(id = c("HR03", "HR04"), geometry = structure(list(
structure(list(list(structure(c(4681404.9717497, 4660789.5055228,
4653707.5146719, 4649417.6856102, 4650218.4304133, 4642887.1080914,
4636001.880903, 4638130.2830531, 4630672.6482304, 4618216.6259566,
4605154.7771129, 4604041.4898082, 4594118.3955698, 4601453.974336,
4625512.2341226, 4632995.7476673, 4630936.0207418, 4634672.4489063,
4642550.4208416, 4643271.7803922, 4667989.5004339, 4677333.5852998,
4698514.5254097, 4706550.7655721, 4730255.159977, 4727301.127819,
4721987.8576822, 4714944.9419401, 4710630.0830272, 4716060.1725087,
4719422.7097735, 4730791.7184246, 4760664.2386483, 4761181.1924784,
4774121.2920662, 4779308.5368569, 4774383.7488244, 4779171.5389403,
4788069.055286, 4799767.4725378, 4802235.1398678, 4799502.7823638,
4807989.5701462, 4809414.3668482, 4813081.7322406, 4810031.5689926,
4817237.2833198, 4815171.9876871, 4817552.444047, 4846574.3460086,
4870047.1348852, 4887573.7485168, 4910786.3738801, 4908879.1782178,
4914796.3717617, 4927371.7873296, 4929499.3861659, 4943730.9291405,
4950718.6813737, 4940446.4383108, 4920595.4028965, 4879438.9383817,
4838212.8736349, 4845048.1363641, 4834560.8316518, 4816355.9591508,
4807977.4769292, 4802984.1968451, 4799132.1108769, 4803857.1680883,
4793328.9122658, 4797545.7578432, 4790387.324717, 4795339.4463412,
4780564.661672, 4768407.3307357, 4746995.6443102, 4729669.638718,
4736099.9440402, 4737387.4547138, 4744143.8223367, 4742293.9670843,
4769269.6338853, 4762671.7474434, 4743127.9197531, 4717672.6244465,
4709482.9634476, 4707387.8919786, 4710229.8207957, 4705993.6051,
4681668.0711335, 4681550.8253001, 4681404.9717497, 2469571.8524468,
2481321.415606, 2458523.7874748, 2449272.8780682, 2438897.1163752,
2438284.5082669, 2422781.4473873, 2418926.8895874, 2414462.7815621,
2437946.1477598, 2451876.6484783, 2465896.8936947, 2492444.7629239,
2491745.7156516, 2487224.7718621, 2490971.3001578, 2496418.20287,
2498159.7068722, 2494429.3004472, 2494333.5855674, 2499017.5900785,
2517639.0231528, 2495368.5990506, 2503501.6452378, 2493771.1676775,
2480763.7475493, 2476651.950761, 2481193.9547111, 2475319.5444323,
2458949.9877284, 2457931.5306729, 2458095.0070284, 2435376.8466673,
2440222.252938, 2442304.3262986, 2432800.850873, 2430045.4784797,
2422025.1304993, 2422657.8559193, 2412903.0969416, 2408110.3193091,
2402794.5017749, 2398224.6047735, 2391577.1648929, 2386104.2008637,
2383099.3384749, 2380363.8037936, 2374363.7382391, 2365756.7861555,
2341238.2275972, 2318986.024962, 2300643.68496, 2290778.8794829,
2283360.2149862, 2272314.0341755, 2260004.8387255, 2258500.0767282,
2251882.761933, 2239655.5626186, 2234753.4147896, 2249058.6706161,
2280366.049377, 2288752.036912, 2292201.2407093, 2293100.1988191,
2282663.0564921, 2285510.2644807, 2285444.7023319, 2299219.0611629,
2299394.7203532, 2308724.3905231, 2308879.7751188, 2315601.3123312,
2309759.8263626, 2312964.4715095, 2323535.336767, 2340392.1183181,
2363926.8629645, 2361734.6731929, 2368030.7347677, 2363291.7705842,
2372862.8783299, 2357030.2024165, 2367280.3689456, 2376330.3981349,
2398007.7669474, 2411137.7946898, 2427411.9337447, 2439546.3651118,
2448439.6717798, 2469358.0996062, 2469630.6323342, 2469571.8524468
), .Dim = c(93L, 2L)), structure(c(4639739.6750726, 4642765.811094,
4639665.9301399, 4639739.6750726, 2440692.0599552, 2438504.6505309,
2444138.6852756, 2440692.0599552), .Dim = c(4L, 2L))), list(
structure(c(4664618.3750015, 4674466.9163483, 4675054.7082127,
4682098.3504266, 4680014.6412814, 4669519.9618955, 4660970.8282448,
4668159.9563802, 4668051.5976083, 4657524.8736319, 4664596.3001291,
4664618.3750015, 2451879.6653522, 2437756.632838, 2412157.9167305,
2402027.1227827, 2399499.2783596, 2408263.7248238, 2434619.0163052,
2431347.488102, 2437766.1868492, 2455007.9464635, 2459754.0870357,
2451879.6653522), .Dim = c(12L, 2L))), list(structure(c(4681404.9717497,
4681668.0711335, 4688448.2211814, 4701638.6577241, 4692877.7268659,
4685560.79219, 4685869.0744863, 4670897.2970111, 4677408.5575198,
4676901.1924166, 4681404.9717497, 2469571.8524468, 2469358.0996062,
2453653.8013503, 2441097.1844415, 2437309.522137, 2441014.8215284,
2447670.7009952, 2449632.430614, 2456680.2670962, 2467759.3352426,
2469571.8524468), .Dim = c(11L, 2L))), list(structure(c(4702920.2744723,
4725652.9674844, 4712128.4539518, 4717239.1438604, 4740982.8946768,
4733702.5011419, 4739222.0941489, 4730779.6749813, 4733170.7052116,
4710696.8794525, 4696480.7246797, 4702920.2744723, 2404000.3011629,
2387695.3114635, 2393696.4815411, 2390106.0734188, 2372294.4781672,
2375411.5731252, 2371045.4605476, 2371670.0920587, 2367943.3447784,
2392439.5205032, 2411114.9411216, 2404000.3011629), .Dim = c(12L,
2L))), list(structure(c(4726698.5350467, 4741471.6282309,
4733519.664328, 4739588.6909075, 4736367.6093496, 4707016.06967,
4709946.7938503, 4726698.5350467, 2335655.2720878, 2324278.926356,
2326563.1197134, 2322039.9508723, 2322772.1693776, 2350755.1956363,
2353407.5676222, 2335655.2720878), .Dim = c(8L, 2L))), list(
structure(c(4881347.5854145, 4878121.2600126, 4841498.1939288,
4843486.1937685, 4881347.5854145, 2270890.3915194, 2265809.1424011,
2267963.3738602, 2276230.1150464, 2270890.3915194), .Dim = c(5L,
2L))), list(structure(c(4907274.0917227, 4863596.0813829,
4839226.8646687, 4862043.3084571, 4868328.5263904, 4907274.0917227,
2252275.7009253, 2247847.7383213, 2253866.0379972, 2258487.4294443,
2252295.6152054, 2252275.7009253), .Dim = c(6L, 2L))), list(
structure(c(5014573.7599384, 5015050.807588, 5023279.372219,
4963657.4237673, 4950932.9886352, 4957005.5664922, 4892115.5770688,
4916054.3231381, 4931608.5403838, 4954856.0811471, 4946488.0750204,
4962248.4460953, 4967064.634373, 4998165.8821143, 5008013.6342012,
5014573.7599384, 2200024.63284, 2192083.0295848, 2185927.9970114,
2220442.7352639, 2224281.5880988, 2216696.3886318, 2242111.7127272,
2239802.630084, 2231244.2030983, 2224658.8064531, 2229853.1835397,
2232980.6019638, 2223168.6861861, 2203437.3540959, 2206028.3808448,
2200024.63284), .Dim = c(16L, 2L))), list(structure(c(4896592.0908832,
4909795.6339469, 4863774.2545106, 4862348.1594577, 4869532.9736551,
4896592.0908832, 2235129.6621158, 2229383.5561162, 2226103.2150051,
2232557.0292488, 2233444.5764789, 2235129.6621158), .Dim = c(6L,
2L)))), class = c("XY", "MULTIPOLYGON", "sfg")), structure(list(
list(structure(c(4735901.1716151, 4741351.1734999, 4764550.6681755,
4761928.980207, 4763299.0547253, 4764748.4752906, 4756287.5678908,
4753400.0163866, 4756763.4914744, 4766659.7722385, 4767522.0004128,
4773915.436675, 4789226.9660126, 4805669.1493577, 4800109.8097679,
4810255.6963994, 4827385.8893674, 4848359.3533164, 4850286.3526126,
4858696.4900945, 4873474.3402934, 4885838.4597189, 4906946.6440309,
4914925.3322593, 4935674.7450365, 4977710.5232465, 4990115.3808876,
4991917.999839, 5004436.5659021, 5009623.0449549, 5007625.9001598,
5013580.7938319, 5008818.7480468, 5018781.6656729, 5019641.784672,
5014675.1111691, 5030998.0535672, 5026320.3224404, 5023243.8778964,
5024713.6482825, 5026974.2280537, 5023139.4347758, 5047125.4042255,
5059718.0933372, 5060947.4131334, 5041089.8801913, 5039312.2849705,
5034795.8358556, 5037785.8503071, 5034268.99946, 5042742.056866,
5030826.9834624, 5033462.499809, 5019939.6619187, 5012646.7541608,
5013996.0922686, 4992124.6831496, 4970013.0938325, 4959986.4708394,
4950513.0337229, 4939160.9256797, 4924125.3744698, 4882603.3106134,
4863942.1431228, 4856677.2149417, 4833894.2445794, 4824332.058776,
4822183.0537321, 4817078.6599674, 4800484.5672889, 4793647.1606918,
4779598.6203156, 4774458.7603405, 4776942.0495754, 4773487.5006867,
4777494.3270629, 4774121.2920662, 4761181.1924784, 4760664.2386483,
4730791.7184246, 4719422.7097735, 4716060.1725087, 4710630.0830272,
4714944.9419401, 4721987.8576822, 4727301.127819, 4730255.159977,
4742233.2265446, 4732879.434036, 4742612.0560869, 4738279.6063048,
4735666.3523642, 4730441.0793027, 4735901.1716151, 2531543.4869999,
2535261.6365574, 2542976.3648409, 2545278.4871249, 2557240.2750943,
2565058.0706851, 2569024.9848123, 2575534.312967, 2583677.4416583,
2583798.4009285, 2589203.2473227, 2592015.892617, 2605342.5752983,
2605574.0857493, 2617566.9916047, 2625262.2117889, 2618351.3261251,
2606260.0956584, 2603109.4583783, 2593197.8787773, 2588532.7445776,
2569452.5404586, 2565907.820542, 2556676.2138283, 2552501.7262843,
2551119.5367503, 2563993.9493555, 2572887.0584876, 2574073.9717903,
2575455.7393524, 2563418.2336765, 2561263.8544271, 2558116.186916,
2556279.9588699, 2547020.2898633, 2537024.7838676, 2532747.3286612,
2528866.5682269, 2528342.0522659, 2522875.2402223, 2520342.9548401,
2515809.6386074, 2504107.3607127, 2504661.947772, 2497458.0581961,
2498871.9000674, 2490525.6360025, 2491695.5083176, 2476998.1898999,
2472765.7390421, 2471166.9996516, 2465486.9613491, 2459084.4515393,
2457365.9738223, 2462095.1960835, 2472080.2476908, 2480502.9942019,
2483453.0791822, 2475736.191394, 2482693.8323636, 2469479.3459692,
2477879.614829, 2476978.7396342, 2487506.8890697, 2476996.4395765,
2479632.8082487, 2466081.9825576, 2453746.3604054, 2452460.0857668,
2463186.2858089, 2475236.4150063, 2474663.5187441, 2467754.1050556,
2462501.4741608, 2456130.6265326, 2446915.4296381, 2442304.3262986,
2440222.252938, 2435376.8466673, 2458095.0070284, 2457931.5306729,
2458949.9877284, 2475319.5444323, 2481193.9547111, 2476651.950761,
2480763.7475493, 2493771.1676775, 2501246.712592, 2513696.834419,
2519238.9170039, 2526272.7357079, 2521664.8325046, 2526256.028918,
2531543.4869999), .Dim = c(94L, 2L)))), class = c("XY",
"MULTIPOLYGON", "sfg"))), class = c("sfc_MULTIPOLYGON", "sfc"
), precision = 0, bbox = structure(c(4594118.3955698, 2185927.9970114,
5060947.4131334, 2625262.2117889), .Names = c("xmin", "ymin",
"xmax", "ymax"), class = "bbox"), crs = structure(list(epsg = 3035L,
proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"), .Names = c("epsg",
"proj4string"), class = "crs"), n_empty = 0L)), .Names = c("id",
"geometry"), row.names = c(NA, -2L), class = c("sf", "data.frame"
), sf_column = "geometry", agr = structure(NA_integer_, class = "factor", .Label = c("constant",
"aggregate", "identity"), .Names = "id"))
NEI = structure(list(id = "BA", geometry = structure(list(structure(list(
list(structure(c(19.0219725000001, 19.171969, 19.3676279990001,
19.303783499, 19.1284175, 19.103181002, 19.1642875000001,
19.3239605020001, 19.3582324980001, 19.623633, 19.5281565000001,
19.2352720000001, 19.5197165, 19.4577080000001, 19.2733945000001,
19.2256605000001, 19.0545335290001, 19.0006863490001, 18.9153873860001,
19.0805207040001, 19.0048281880001, 18.9102436480001, 18.6884211570001,
18.6631543020001, 18.50043078, 18.4599684510001, 18.565898029,
18.5237265760001, 18.4372755000001, 18.417495023, 18.3663980010001,
18.342735746, 18.2440605, 18.2407502080001, 18.1572590500001,
17.957566465, 17.891265999, 17.8090535000001, 17.7624947710001,
17.7432932550001, 17.6491623530001, 17.6491565000001, 17.6161355650001,
17.6113094340001, 17.594219876, 17.5813350000001, 17.7120725000001,
17.6412180010001, 17.4508455, 17.2654190000001, 17.2824995000001,
17.0079865, 16.544348, 16.1724470000001, 16.2194785010001,
16.2187025000001, 16.1314050000001, 16.1724394990001, 16.1319370000001,
16.1205805010001, 16.0186899990001, 16.031724998, 15.8940600000001,
15.834873, 15.781581, 15.7367220010001, 15.7858065000001,
15.764635998, 15.8362885000001, 16.0149785, 16.0897560000001,
16.2927725, 16.3545210000001, 16.3943880000001, 16.529669,
16.8154390000001, 16.9338705, 16.9383715, 17.1435915, 17.2456665000001,
17.271448001, 17.4827439990001, 17.6697865, 17.8656290010001,
17.9941825020001, 18.1200030000001, 18.1966040000001, 18.2568340000001,
18.4777605, 18.5319975, 18.6677210000001, 18.8497695000001,
19.0219725000001, 44.855506, 44.920824001, 44.8785735, 44.689073,
44.519913, 44.368331, 44.285345, 44.271836999, 44.185297,
44.050928, 43.9567225, 44.0110395, 43.7101305, 43.5979170000001,
43.5970585, 43.5284235, 43.501701146, 43.558701364, 43.503957057,
43.314686475, 43.253429821, 43.360430057, 43.254683833, 43.036829651,
42.994723683, 42.821522782, 42.662171666, 42.584645607, 42.559768499,
42.5755013420001, 42.616142499, 42.613542384, 42.6026995,
42.604642478, 42.653647677, 42.770857391, 42.8097725, 42.917099,
42.908894872, 42.905511367, 42.888924786, 42.8889235, 42.913023819,
42.916546169, 42.929018978, 42.938423, 42.9730035, 43.089164501,
43.1768685, 43.372799, 43.4677620000001, 43.5762404990001,
43.9741785000001, 44.2092095, 44.223788999, 44.348137, 44.3780675,
44.4028015, 44.4546815, 44.5154455, 44.562702, 44.653316501,
44.7493515, 44.713775501, 44.7499275, 44.935631, 44.9941445,
45.1640165, 45.222477, 45.217575, 45.1043815, 44.9989320000001,
45.0033450000001, 45.112419, 45.226574, 45.184578, 45.276321499,
45.2273715, 45.162589501, 45.1460075000001, 45.1889875000001,
45.1101495, 45.1335680000001, 45.0453225, 45.148647499, 45.0802305,
45.0785905, 45.139408, 45.062134, 45.0906295, 45.061779,
44.8545455, 44.855506), .Dim = c(93L, 2L)))), class = c("XY",
"MULTIPOLYGON", "sfg"))), class = c("sfc_MULTIPOLYGON", "sfc"
), precision = 0, bbox = structure(c(15.7367220010001, 42.559768499,
19.623633, 45.276321499), .Names = c("xmin", "ymin", "xmax",
"ymax"), class = "bbox"), crs = structure(list(epsg = NA_integer_,
proj4string = "+proj=longlat +ellps=GRS80 +no_defs"), .Names = c("epsg",
"proj4string"), class = "crs"), n_empty = 0L)), .Names = c("id",
"geometry"), row.names = c(NA, -1L), class = c("sf", "data.frame"
), sf_column = "geometry", agr = structure(NA_integer_, class = "factor", .Label = c("constant",
"aggregate", "identity"), .Names = "id"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.