Map_making_BHA

knitr::opts_chunk$set(echo = TRUE)

### Tell R where your data is stored

setwd("C:/Users/N11094427/OneDrive - Queensland University of Technology/Documents/tmp/BHA_map") # substitude to your own path

library(ggplot2)
library(ozmaps)
library(raster)

## Loading required package: sp

library(gridExtra)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sf)

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# We need to define the area of interest:

BHA_property <- st_read("qld_remnant_re_tierV1.shp")

## Reading layer `qld_remnant_re_tierV1' from data source 
##   `C:\Users\N11094427\OneDrive - Queensland University of Technology\Documents\tmp\BHA_map\qld_remnant_re_tierV1.shp' 
##   using driver `ESRI Shapefile'

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: C:\Users\N11094427\OneDrive - Queensland University of
## Technology\Documents\tmp\BHA_map\qld_remnant_re_tierV1.shp contains polygon(s)
## with rings with invalid winding order. Autocorrecting them, but that shapefile
## should be corrected using ogr2ogr for example.

## Simple feature collection with 1015 features and 16 fields (with 5 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 141.5358 ymin: -26.04321 xmax: 149.8036 ymax: -20.80835
## Geodetic CRS:  GCS_GDA2020

unique(BHA_property$property)

## [1] "Pony Hill"     "Alpha station" "Pegunny"       "Oakwood"      
## [5] "Coolreagh"     "Yorkshire"     NA

# Filter the properties of interest
Pony_hill <- filter(BHA_property, property == "Pony Hill")
Alpha_station <- filter(BHA_property, property == "Alpha Station")
Pegunny <- filter(BHA_property, property == "Pegunny")
Coolreagh <- filter(BHA_property, property == "Coolreagh")
Oakwood <- filter(BHA_property, property == "Oakwood")

library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(ggspatial)
library(prismatic) #Approximation to hexcode colors in console

# remotes::install_github('qdread/qdrmapbox') #this is a package for using the free statellite imageries offered by Mapbox 

library(qdrmapbox)

set_mapbox_api_key('C:/Users/N11094427/OneDrive - Queensland University of Technology/Documents/tmp/BHA_map/mapboxapikey.txt') # you need to register a free account and get a token to download free imageries from Mapbox

## Warning in readLines(keyfile): incomplete final line found on
## 'C:/Users/N11094427/OneDrive - Queensland University of
## Technology/Documents/tmp/BHA_map/mapboxapikey.txt'

download_dir <- 'C:/Users/N11094427/OneDrive - Queensland University of Technology/Documents/tmp/BHA_map/mdtiles'  # this is the place you would like to store your downloaded satellite imageries


zoom <- 6  # zoom is equal to 6 means the whole earth is covered by 2^6 = 64 tiles (resolution ~ 2.5 km), this a numbering system used by common map servers provider like Openstreetmap and google map

library(ozmaps)

# Load the Queensland map
QLD <- ozmap_states[ozmap_states$NAME == "Queensland",]

# Transform the CRS of Queensland map to match the property data
QLD <- st_transform(QLD, crs(BHA_property))

#Basemap: QLD
QLDExtent <- st_bbox(QLD) # we need the extend of interested area

upper_left <- c(QLDExtent[["ymax"]], QLDExtent[["xmin"]]) # the upper left of the extent of QLD
lower_right <- c(QLDExtent[["ymin"]], QLDExtent[["xmax"]]) # the lower right of the extent of QLD

QLD_md_tile_index <- find_tile_numbers(zoom = zoom, upper_left = upper_left, lower_right = lower_right)

QLD_md_tile_df <- download_mapbox_tiles(tile_numbers_mat = QLD_md_tile_index, download_dir = download_dir, resolution = 'high', jpg_quality = 90) # download satellite imageries, resolution choice: "low" and "high", jpg_quality choice: "70","80","90"

georeference_all_tiles(QLD_md_tile_df)

build_virtual_raster(QLD_md_tile_df, file.path(download_dir, 'QLD_mdimage.vrt'))

QLD_mdimage_raster <- stack(file.path(download_dir, 'QLD_mdimage.vrt'))

QLD_df <- as.data.frame(QLD_mdimage_raster, xy= TRUE)

library(dplyr)

QLD_df  <- QLD_df  %>% rename(Red = QLD_mdimage_1,   #Rename bands
                    Green = QLD_mdimage_2,
                    Blue = QLD_mdimage_3)

# Pny hill map

Pony_hill_Extent <- st_bbox(Pony_hill)

Pony_hill_upper_left <- c(Pony_hill_Extent[["ymax"]]+ 0.1, Pony_hill_Extent[["xmin"]]- 0.1) # the upper left of the extent of QLD
Pony_hill_lower_right <- c(Pony_hill_Extent[["ymin"]]- 0.1, Pony_hill_Extent[["xmax"]]+ 0.1) # the lower right of the extent of QLD

QLD_map <-  ggplot()+                 
  geom_raster(data = QLD_df , aes(x = x, y =y),
              fill = rgb(r = QLD_df$Red,
                         g = QLD_df$Green,
                         b = QLD_df$Blue,
                         maxColorValue = 255),
                         show.legend = FALSE) +
  xlim(QLDExtent[["xmin"]]-0.5,QLDExtent[["xmax"]]+0.5)+
  ylim(QLDExtent[["ymin"]]-0.5,QLDExtent[["ymax"]]+0.5)+
  geom_sf(data = st_geometry(QLD), color = 'white', fill = NA, size = 2)+
  geom_sf(data = Pony_hill, color = "white", fill = NA) +
  annotate("rect", 
           xmin =Pony_hill_Extent["xmin"]-0.2, xmax = Pony_hill_Extent["xmax"]+0.2, 
           ymin = Pony_hill_Extent["ymin"]-0.2, ymax = Pony_hill_Extent["ymax"]+0.2, 
           color = "yellow", fill = NA, size = 0.5) +  #the yellow rectangle around Pony hill in the QLD map
  theme(axis.title.x=element_blank(),axis.title.y =element_blank(),
        axis.text.x=element_blank(), axis.text.y=element_blank(),
        axis.ticks=element_blank(),plot.margin = unit(c(0, 0, 0, 0), "cm"))+
  theme_void()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

QLD_map

## Warning: Removed 2178170 rows containing missing values or values outside the scale
## range (`geom_raster()`).

PH_md_tile_index <- find_tile_numbers(zoom = 11, upper_left = Pony_hill_upper_left, lower_right = Pony_hill_lower_right) # resolution ~ 80m

PH_md_tile_index <- download_mapbox_tiles(tile_numbers_mat = PH_md_tile_index, download_dir = download_dir, resolution = 'high', jpg_quality = 90) # download satellite imageries, resolution choice: "low" and "high", jpg_quality choice: "70","80","90"

georeference_all_tiles(PH_md_tile_index)

build_virtual_raster(PH_md_tile_index, file.path(download_dir, 'PH_high_resolution_mdimage.vrt'))

PH_high_mdimage_raster <- stack(file.path(download_dir, 'PH_high_resolution_mdimage.vrt'))

PH_df <- as.data.frame(PH_high_mdimage_raster, xy= TRUE)

library(dplyr)

PH_df <- PH_df %>% rename(Red = PH_high_resolution_mdimage_1,   #Rename bands
                    Green = PH_high_resolution_mdimage_2,
                    Blue = PH_high_resolution_mdimage_3)

Pony_hill_map <- 
  ggplot()+                   #plot map
  geom_raster(data = PH_df , aes(x = x, y =y),
              fill = rgb(r = PH_df$Red,
                         g = PH_df$Green,
                         b = PH_df$Blue,
                         maxColorValue = 255),
                         show.legend = FALSE) +
  xlim(Pony_hill_Extent[["xmin"]]-0.2,Pony_hill_Extent[["xmax"]]+0.05)+
  ylim(Pony_hill_Extent[["ymin"]]-0.1,Pony_hill_Extent[["ymax"]]+0.05)+
  geom_sf(data = st_geometry(Pony_hill),color = 'white', fill = NA, size =0.5)+
  ggtitle("Pony Hill") +
  annotation_scale(location = "br", text_col= 'white',
                   pad_x = unit(1, "cm"),, pad_y = unit(1, "cm"),)+
  annotation_north_arrow(location = "tl", which_north = "true",
                         height = unit(1, "cm"), width = unit(1, "cm"),
                         pad_x = unit(1, "cm"),, pad_y = unit(1.2, "cm"),
                         style = north_arrow_fancy_orienteering(text_col = 'white',
                                                                line_col = 'white',
                                                                fill = 'white'))+
  annotate("text", x = 149.255, y = -26.142, label = '\u00a9 Mapbox \u00a9 OpenStreetMap', 
           hjust = 1, vjust = -1, color = 'white', size = 2) +
  theme_bw()+
  theme(axis.title.x=element_blank(),axis.title.y =element_blank(),
        plot.margin = unit(c(1, 1, 1, 1), "cm"))

Pony_hill_map

## Warning: Removed 1859160 rows containing missing values or values outside the scale
## range (`geom_raster()`).

library(cowplot)

# Combine the main map and the inset map
combined_map <- ggdraw() +
  draw_plot(Pony_hill_map) +
  draw_plot(QLD_map, x = 0.18, y = 0.17, width = 0.3, height = 0.3)

## Warning: Removed 1859160 rows containing missing values or values outside the scale
## range (`geom_raster()`).

## Warning: Removed 2178170 rows containing missing values or values outside the scale
## range (`geom_raster()`).

# Print the combined map
print(combined_map)

ggsave("Pony_hill.png", width = 5000, height = 3500, units = "px", dpi = 800)