How to extract also a zip code with Overpass turbo?

Question

I use this code in Overpass turbo:

[out:csv(comune, ISTAT, name, highway, ::lat,::lon;false;";")][timeout:600];
//provincia da cui estrarre i dati
area[boundary=administrative]["admin_level"=6][name="Napoli"]->.searchArea;
relation[boundary=administrative]["admin_level"=8](area.searchArea);
foreach (
  map_to_area->.comune;
  make stat comune=comune.set(t["name"]),ISTAT=comune.set(t["ref:ISTAT"]);
  out;
  way[highway~"residential|unclassified|tertiary|secondary|primary"][name](area.comune);
  out center;
);

Is it possible to extract also the zip postal code?

Answer 1

There is no much municipalities with postal code in your area of interest:

            name postal_code
40663       Meta       80062
40730     Pompei       80045
40876 Casavatore       80020

My attempt in R to combine OSM with Wikidata. For OSM retrieval will use osmdata package. As a first step we are looking for Napoli bounding box, then for administrative boundary on level 6:

napoli_bb <- osmdata::getbb("Napoli")

napoli <- osmdata::opq(napoli_bb, timeout = 60*20) |>
  osmdata::add_osm_feature(key = "boundary", value = "administrative") |>
  osmdata::add_osm_feature(key = "admin_level", value = "6") |>
  osmdata::add_osm_feature(key = "name", value = "Napoli") |>
  osmdata::osmdata_sf() |>
  osmdata::unique_osmdata() |>
  osmdata::unname_osmdata_sf()

napoli_boundary <- napoli$osm_multipolygons |>
  subset(name == "Napoli" & boundary == "administrative" & admin_level == 6)

Second step, we are using the boundary bounding box to download all administrative boundaries and highways

napoli <- osmdata::opq(bbox = sf::st_bbox(napoli_boundary), timeout = 60*20) |>
  osmdata::add_osm_features(features = 
                              c("\"boundary\" = \"administrative\"",
                                "\"highway\"")) |>
  osmdata::osmdata_sf() |>
  osmdata::unique_osmdata() |>
  osmdata::unname_osmdata_sf()

Then just filter out those municipalities, which are within Napoli boundaries.

municipalities <- napoli$osm_multipolygons |>
  subset(boundary == "administrative" & admin_level == 8) |>
  sf::st_filter(napoli_boundary) |>
  subset(select = c("name", "ref:ISTAT", "wikidata"))

With help of WikidataR package and get_postal_code_from_wikidata() function we are querying wikidata and retrieve the postal code(s). Some municipalities have more than one. It's added to municipalities data frame.

get_postal_code_from_wikidata <- function(q = "") {
  if(q != "" & !is.na(q)) {
    w <- WikidataR::get_item(q)
    x <- WikidataR::extract_claims(w, "P281")
    a <- x[[1]]$P281.x$mainsnak$datavalue$value
    if(length(a > 1)) {
      a <- a |>
        paste0(str = ",", collapse = " ")
    }
  } else {
    a <- ""
  }  
  return(stringr::str_replace(a, ",$", ""))
}

municipalities <- municipalities |>
  dplyr::rowwise() |>
  dplyr::mutate(postal_code = get_postal_code_from_wikidata(wikidata))

That's the result:

municipalities |>
  sf::st_drop_geometry()
#> # A tibble: 134 × 4
#> # Rowwise: 
#>    name              `ref:ISTAT` wikidata postal_code        
#>  * <chr>             <chr>       <chr>    <chr>              
#>  1 Anacapri          063004      Q71617   80071              
#>  2 Capri             063014      Q71902   80073              
#>  3 Sorrento          063080      Q72672   80067              
#>  4 Sant'Agnello      063071      Q72565   80065              
#>  5 Furore            065053      Q80948   84010              
#>  6 Piano di Sorrento 063053      Q72346   80063              
#>  7 Praiano           065102      Q81355   84010              
#>  8 Vico Equense      063086      Q72732   80069, 80060, 80066
#>  9 Amalfi            065006      Q80563   84011              
#> 10 Positano          065100      Q81345   84017              
#> # ℹ 124 more rows

Now, as an example, filtering out the roads from our downloaded data, for one municipality: Barano d'Ischia. Please note, the roads are grouped (their geometries are 'unioned') by name and type of highway. center_lon and _lat are coordinates of their centroids.

municipalities[23,,]$name
#> [1] "Barano d'Ischia"


roads <- napoli$osm_lines |>
  subset(highway %in% c("residential", "unclassified", "tertiary", 
                        "secondary", "primary"), select = c(name, highway)) |>
  sf::st_filter(municipalities[23, ,]) |>
  dplyr::group_by(name, highway) |>
  dplyr::summarise(.groups = "keep") |>
  dplyr::mutate(center_lon = sf::st_coordinates(sf::st_centroid(geometry))[[1]],
                center_lat = sf::st_coordinates(sf::st_centroid(geometry))[[2]])

roads |>
  sf::st_drop_geometry()
#> # A tibble: 22 × 4
#> # Groups:   name, highway [22]
#>    name                         highway    center_lon  center_lat
#>  * <chr>                        <chr>           <dbl>       <dbl>
#>  1 Piazza San Rocco             secondary        13.9        40.7
#>  2 Via Angelo Migliaccio        secondary        13.9        40.7
#>  3 Via Astiere                  unclassified     13.9        40.7
#>  4 Via Fasolara                 residential      13.9        40.7
#>  5 Via Nuova dei Conti          tertiary         13.9        40.7
#>  6 Via Nuova dei Conti          unclassified     13.9        40.7
#>  7 Via Pendio del Gelso         tertiary         13.9        40.7
#>  8 Via Piano                    residential      13.9        40.7
#>  9 Via Provinciale Maronti      tertiary         13.9        40.7
#> 10 Via Provinciale Regina Elena unclassified     13.9        40.7
#> # ℹ 12 more rows

And small map with our data:

tmap::tm_shape(napoli_boundary) +
  tmap::tm_basemap("OpenStreetMap") +
  tmap::tm_borders(lwd = 3) +
  tmap::tm_shape(municipalities) +
  tmap::tm_polygons(col = "grey95") +
  tmap::tm_shape(roads) +
  tmap::tm_lines()

^{Created on 2024-01-23 with reprex v2.1.0}

You can ran it for each municipality with for() loop or make a function and bind the rows. Or use kind of municipalities |> sf::st_join(napoli$osm_lines) and group by municipality, but it will not cut the highways on the boundaries.

Please note that coordinates are given with few more digits after point, it's just a matter of output format:

roads[1, ]$center_lon
#> [1] 13.91921