Getting and cleaning data, example of Chinese airports - Part 4/5

Reorganising columns

We update column order to get a clearer view when comparing data.

# reorder
IATA <- dplyr::select(IATA,
                      IATA.FAA, ICAO.x, ICAO.y,
                      name.x, name.y, city.x, city.y,
                      latitude.deg.x, latitude.deg.y,
                      longitude.deg.x, longitude.deg.y,
                      elevation.m.x, elevation.m.y,
                      category.x, category.y,
                      region.name, type.airport, scheduled.service, airport.wikipedia.link
) 

Fixing ICAO codes

Let’s clean the ICAO codes, we will use the one from OurAiports. All missing values of OpenFlight are included in Ouraiports. Missing values of OurAirports are not found in Openflight data set (see previous missingness map)

IATA$ICAO <- IATA$ICAO.x
IATA <- dplyr::select(IATA, -ICAO.x, -ICAO.y)
missmap(IATA)

Fixing city name

No missing data in city.y, we replace missing city.x by associated city.y. There seems to be small typos in both data set (always for Chinese cities), we compare city names by calculating the distance between the character vectors. We keep city.x which seems better.

IATA[is.na(IATA$city.x),]$city.x <- IATA[is.na(IATA$city.x),]$city.y 

dist <- stringdist(IATA$city.x, IATA$city.y)
hist(dist, breaks = 40)

kable(dplyr::select(IATA[dist > 0,], city.x, city.y, name.x))

IATA$city <- IATA$city.x
IATA <- dplyr::select(IATA, -city.x, -city.y)
missmap(IATA)

Fixing airport names

No missing data in name.y or name.x, we need to compare both by string distance. name.x is more detailed and better written. Except for Taoxian Airport.

dist <- stringdist(IATA$name.x, IATA$name.y)
hist(dist, breaks = 40)

kable(head(dplyr::select(IATA[dist > 0,], name.x, name.y)))

IATA[IATA$name.x == "Taoxian Airport",]$name.x <- "Shenyang Taoxian International Airport"
IATA$name <- IATA$name.x
IATA <- dplyr::select(IATA, -name.x, -name.y)
missmap(IATA)

Fixing longitude.deg, latitude.deg

Values for longitude and latitude are different in the two data sets, we study in details the distance over 0.01 degree which represent roughly 500 meters.

dist1 <- abs(IATA$longitude.deg.x - IATA$longitude.deg.y)
dist2 <- abs(IATA$latitude.deg.x - IATA$latitude.deg.y)

par(mfrow = c(1,2))
hist(log(dist1), breaks = 40)
hist(log(dist2), breaks = 40)

toCheck <- dplyr::select(IATA[dist1 > 0.05 & dist2 > 0.05,],
                         IATA.FAA, ICAO, airport.wikipedia.link,
                         latitude.deg.x, longitude.deg.x,
                         latitude.deg.y,  longitude.deg.y)

kable(toCheck)

Let’s now plot these airports.

g <- ggplot() + theme_bw()
## airports with elevation
g <- g + geom_point(data = toCheck, aes(x = longitude.deg.x, y = latitude.deg.x), col = "firebrick", pch = 3)
## airport with 0 m  as elevation in red
g <- g + geom_point(data = toCheck, aes(x = longitude.deg.y, y = latitude.deg.y), col = "blue", pch = 1)
g <- g + geom_text(data = toCheck, aes(x = longitude.deg.x, y = latitude.deg.x, label=IATA.FAA), size = 3, vjust = -1)
print(g)

The map above shows there are significant error between the two data sets, now let’s sample one case to check which data set is probably the most accurate. We used for that the OpenStreetMap package to display side by side satellite images around the latitude and longitude. This worked but for some unknown reasons some tiles cannot be downloaded from travis. So I changed strategy, using the API of mapbox instead but with longer and more complicated code.

# set-up size of the bounding box to be visualized.
delta_lon <- 0.01
delta_lat <- 0.005
zoom_level <- 15


mapAirport <- function(
  my_lon, my_lat,
  my_delta_lat = 0.01,
  my_delta_lon = 0.005,
  my_zoom_level = 14,
  airport_id = 1,
  my_api_mapbox
) {
  ## Define bounding box
  myBbox <- c(xmin = my_lon - my_delta_lon,
              ymin = my_lat - my_delta_lat,
              xmax = my_lon + my_delta_lon,
              ymax = my_lat + my_delta_lat
  )
  # make the gridding
  myBboxTile <- slippymath::bbox_tile_query(bbox = myBbox)
  myGrid <- slippymath::bbox_to_tile_grid(myBbox, zoom = my_zoom_level)
  
  
  ## Download satellite data
  setwd(dir = "./sat/")
  
  mapbox_query_string <-
    paste0(
      "https://api.mapbox.com/v4/mapbox.satellite/{zoom}/{x}/{y}.jpg90",
      #"https://api.mapbox.com/v4/mapbox.satellite/{zoom}/{x}/{y}@2x.jpg90",
      "?access_token=",
      my_api_mapbox
    )
  
  
  
  ## Function to download Tiles
  myTiles <-
    pmap(.l = myGrid$tiles,
         zoom = myGrid$zoom,
         
         .f = function(x, y, zoom){
           outfile <- glue("{x}_{y}.jpg")
           curl_download(url = glue(mapbox_query_string),
                         destfile = outfile)
           outfile
         }
    )
  
  
  mySat <- slippymath::compose_tile_grid(myGrid, myTiles)
  #raster::plot(mySat)
  
  ## Go back to root directory
  setwd(dir = "../")
  #raster_to_png(mySat, "../../static/img/md/test_mapbox.png")
  raster_to_png(mySat, paste0("./img/airport_",airport_id,".png"))
  #knitr::include_graphics("/img/md/test_mapbox.png")
  #knitr::include_graphics(paste0("/blog/img/airport_",airport_id,".png"))
  return(paste0("/blog/img/airport_",airport_id,".png"))
  #overlay_image2 <- png::readPNG("final.png")
}

list_img <- NULL

for( i in 1:nrow(toCheck)){
  cat(paste0("<br/><br/><b>",toCheck[i,]$IATA.FAA,
             " - <a href='",toCheck[i,]$airport.wikipedia.link,"'>",
             toCheck[i,]$airport.wikipedia.link,"</a></b></br>\n"))
  a <- mapAirport(
    my_lon = toCheck[i,]$longitude.deg.x,
    my_lat = toCheck[i,]$latitude.deg.x,
    my_delta_lat = delta_lat,
    my_delta_lon = delta_lon,
    my_zoom_level = zoom_level,
    airport_id = 2*i,
    my_api_mapbox = api_mapbox
  )
  #cat(paste0("![",toCheck[i,]$IATA.FAA,"_OA](", a, " =300x)\n"))
  cat(paste0("<img src=",a," alt='",toCheck[i,]$IATA.FAA,"_OA'"," width='300'/>\n"))
  b <- mapAirport(
    my_lon = toCheck[i,]$longitude.deg.y,
    my_lat = toCheck[i,]$latitude.deg.y,
    my_delta_lat = delta_lat,
    my_delta_lon = delta_lon,
    my_zoom_level = zoom_level,
    airport_id = 2*i+1,
    my_api_mapbox = api_mapbox
  )
  #cat(paste0("![",toCheck[i,]$IATA.FAA,"_OF](", b, " =300x)\n"))
  cat(paste0("<img src=",b," alt='",toCheck[i,]$IATA.FAA,"_OF'"," width='300'/>\n"))
}

CGQ - https://en.wikipedia.org/wiki/Changchun_Longjia_International_Airport


DAX - https://en.wikipedia.org/wiki/Daxian_Airport


DNH - https://en.wikipedia.org/wiki/Dunhuang_Airport


DOY - https://en.wikipedia.org/wiki/Dongying_Shengli_Airport


DSN - https://en.wikipedia.org/wiki/Ordos_Ejin_Horo_Airport


FUO - https://en.wikipedia.org/wiki/Foshan_Shadi_Airport


GOQ - https://en.wikipedia.org/wiki/Golmud_Airport


GXH - https://en.wikipedia.org/wiki/Gannan_Xiahe_Airport


HZG - https://en.wikipedia.org/wiki/Hanzhong_Chenggu_Airport


HZH - https://en.wikipedia.org/wiki/Liping_Airport


INC - https://en.wikipedia.org/wiki/Yinchuan_Hedong_International_Airport


IQM - https://en.wikipedia.org/wiki/Qiemo_Yudu_Airport


JIU - https://en.wikipedia.org/wiki/Jiujiang_Lushan_Airport


JNG - https://en.wikipedia.org/wiki/Jining_Qufu_Airport


KMG - https://en.wikipedia.org/wiki/Kunming_Changshui_International_Airport


KRY - https://en.wikipedia.org/wiki/Karamay_Airport


LYA - https://en.wikipedia.org/wiki/Luoyang_Beijiao_Airport


SHE - https://en.wikipedia.org/wiki/Shenyang_Taoxian_International_Airport


SWA - https://en.wikipedia.org/wiki/Jieyang_Chaoshan_International_Airport


THQ - https://en.wikipedia.org/wiki/Tianshui_Maijishan_Airport


TLQ - https://en.wikipedia.org/wiki/Turpan_Jiaohe_Airport


UYN - https://en.wikipedia.org/wiki/Yulin_Yuyang_Airport


WNH - https://en.wikipedia.org/wiki/Wenshan_Puzhehei_Airport


XUZ - https://en.wikipedia.org/wiki/Xuzhou_Guanyin_Airport


YNT - https://en.wikipedia.org/wiki/Yantai_Penglai_International_Airport

As a conclusion for this part, OA seems to be more accurate than OF regarding the airport longitude and latitude. We decide to keep OA data and drop OF data. note we should confirm that openstreetmap satellite is accurate in China and there is not shift phenomenon.

IATA$latitude.deg <- IATA$latitude.deg.x
IATA$longitude.deg <- IATA$longitude.deg.x
IATA <- dplyr::select(IATA, -latitude.deg.x, -latitude.deg.y, -longitude.deg.x, -longitude.deg.y)

Nevertheless as shown in the case of ZLXH airport the data of OA seems inaccurate. Further check is probably needed. note 2 Comparison is done here manually, we will see in another post how to make this detection using machine learning algorithms.

Fixing categories

They are same in both data sets except the following. We decide to keep OF categories.

## mismatch in category
IATA[!(IATA$category.x == IATA$category.y),]$name

## [1] "Altay Air Base"    "Shigatse Air Base"

IATA$category <- IATA$category.y
IATA <- dplyr::select(IATA, -category.x, -category.y)
missmap(IATA)

Fixing elevation.m data

As we have decided to drop OF latitude and longitude we should also drop the elevation from OF data set.

IATA$elevation.m <- IATA$elevation.m.x
IATA <- dplyr::select(IATA, -elevation.m.y, -elevation.m.x)

We have two actions to check and fix elevation data : check if the elevation is consistent with latitude and longitude and find missing elevation data.

We can use for this purpose the jawq API which returns in JSON format the elevation value for a given set of latitude / longitude.

# function which return evelevation based on longitude / latitude
# https://developers.google.com/maps/documentation/elevation/intro

getElevation <- function(lat,long) {
  lat <- as.numeric(lat)
  long <- as.numeric(long)
  # elevation API
  #baseUrl <- "https://api.jawg.io/elevations?locations="
  baseUrl <- "https://maps.googleapis.com/maps/api/elevation/json?locations="

  Url <- paste(baseUrl, lat, "," , long, "&key=", api_elevation, sep = "")
  
  jsonData <- fromJSON(Url)
  #Sys.sleep(1)
  elevation <- jsonData$results$elevation
  elevation
  #Url
}
## get elevation for all aiports
#IATA$elevationCheck <- 1
elevationCheck <- apply(IATA, 1, function(x) getElevation(x[9],x[10]))

IATA$elevationCheck <- elevationCheck

We are now going to check the gap between Google data and OA data set.

## replace missing elevation by google values
IATA[is.na(IATA$elevation.m),]$elevation.m <- IATA[is.na(IATA$elevation.m),]$elevationCheck
dist <- (abs(IATA$elevation.m - IATA$elevationCheck))
hist(dist, breaks = 40)

summary(dist)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0000    0.1071    2.1109   50.4120    6.2938 2648.6264

dplyr::select(IATA[dist > 30,], IATA.FAA, ICAO, elevation.m, elevationCheck, longitude.deg, latitude.deg)

##     IATA.FAA ICAO elevation.m elevationCheck longitude.deg latitude.deg
## 26       DCY ZUDC    4419.749      4389.1841     100.05333     29.32306
## 46       HEK ZYHE    2605.062       316.7160     127.30888     50.17162
## 79       KCA ZWKC    1076.230      1031.4375      82.87292     41.67786
## 116      NNY ZHNY     256.536       123.4809     112.61500     32.98080
## 122      PZI ZUZH     494.748      1963.4302     101.79852     26.54000
## 123      RKZ ZURK    1155.023      3803.6492      89.31140     29.35190
## 137      TGO ZBTL     731.433       183.6307     122.20000     43.55670
## 178      ZHY ZLZW    2504.891      1237.2886     105.15445     37.57312

We have 8 airports with more than 30 meters of elevation error. After checking, we find out that Jawg data is more accurate and use it for these airports.

IATA[dist > 30,]$elevation.m <- IATA[dist > 30,]$elevationCheck
IATA <- dplyr::select(IATA, -elevationCheck)
missmap(IATA)

We use information from [http://www.geoplaner.com/#Coordinate-Converter] for setting

Source code

The source code of this post is available on github

Session information

sessionInfo()

## R version 3.6.1 (2017-01-27)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /home/travis/R-bin/lib/R/lib/libRblas.so
## LAPACK: /home/travis/R-bin/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
##  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
##  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
## [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] curl_4.3            purrr_0.3.3         png_0.1-7          
##  [4] raster_3.0-7        glue_1.3.1          slippymath_0.3.1   
##  [7] knitr_1.26          stringdist_0.9.5.5  OpenStreetMap_0.3.4
## [10] Amelia_1.7.6        Rcpp_1.0.3          jsonlite_1.6       
## [13] ggplot2_3.2.1       dplyr_0.8.3         reshape2_1.4.3     
## [16] maptools_0.9-9      sp_1.3-2           
## 
## loaded via a namespace (and not attached):
##  [1] highr_0.8        compiler_3.6.1   pillar_1.4.2     plyr_1.8.5      
##  [5] tools_3.6.1      digest_0.6.23    evaluate_0.14    tibble_2.1.3    
##  [9] lifecycle_0.1.0  gtable_0.3.0     lattice_0.20-38  pkgconfig_2.0.3 
## [13] rlang_0.4.2      parallel_3.6.1   rgdal_1.4-8      yaml_2.2.0      
## [17] blogdown_0.17.1  xfun_0.11        rJava_0.9-11     withr_2.1.2     
## [21] stringr_1.4.0    tidyselect_0.2.5 R6_2.4.1         foreign_0.8-72  
## [25] rmarkdown_1.18   bookdown_0.16    farver_2.0.1     magrittr_1.5    
## [29] codetools_0.2-16 scales_1.1.0     htmltools_0.4.0  assertthat_0.2.1
## [33] colorspace_1.4-1 labeling_0.3     stringi_1.4.3    lazyeval_0.2.2  
## [37] munsell_0.5.0    crayon_1.3.4