Shapefiles and Raster of Relief do not match for Hawaii

gistfile1.txt

### This method works for Switzerland and The USA (shown here, at the top) but not Hawaii (shown at the bottom)
library(sf)
library(raster)
library(dplyr)
library(ggplot2)

## USA
# download shapefile: https://catalog.data.gov/dataset/tiger-line-shapefile-2017-nation-u-s-current-state-and-equivalent-national 
usa_geo <- read_sf("data/shapefiles/usa/tl_2017_us_state.shp")

# download relief: http://www.naturalearthdata.com/downloads/50m-raster-data/50m-shaded-relief/
relief <- raster("data/relief/world/SR_50M/SR_50M.tif") %>%
  as("SpatialPixelsDataFrame") %>%
  as.data.frame() %>%
  rename(., value = "SR_50M") %>% 
  filter(value != 255)

ggplot(
    # define main data source
    data = usa_geo
) +
    # first: draw the relief
    geom_raster(
        data = relief,
        inherit.aes = FALSE,
        aes(
            x = x,
            y = y,
            alpha = value
        )
    ) +
    # use the "alpha hack" (as the "fill" aesthetic is already taken)
    scale_alpha(name = "",
                range = c(0.4, 0.01),
                guide = F) + # suppress legend
    # add main fill aesthetic
    # use thin white stroke for municipality borders
    geom_sf(color = "red",
            size = 0.1
    )
    
## Hawaii
# https://earthworks.stanford.edu/catalog/stanford-qh711pf3383
relief <- raster("data/relief/100m/srgrhii0100a.tif") %>%
  as("SpatialPixelsDataFrame") %>%
  as.data.frame() %>%
  rename(., value = "srgrhii0100a") %>% 
  filter(value != 255)

# http://geoportal.hawaii.gov/datasets/045b1d5147634e2380566668e04094c6_3
hawaii_geo <- read_sf("data/shapefiles/coastline/Coastline.shp")

ggplot(
    # define main data source
    data = hawaii_geo
) +
    # first: draw the relief
    geom_raster(
        data = relief,
        inherit.aes = FALSE,
        aes(
            x = x,
            y = y,
            alpha = value
        )
    ) +
    # use the "alpha hack" (as the "fill" aesthetic is already taken)
    scale_alpha(name = "",
                range = c(0.4, 0.01),
                guide = F) + # suppress legend
    # add main fill aesthetic
    # use thin white stroke for municipality borders
    geom_sf(color = "red",
            size = 0.1
    )

Hi @moldach, can you share all of the required files in a one archive (e.g., a zip file)? It would help me to reproduce your problem without going to every single webpage and downloading data one by one...

Hi @moldach, can you share all of the required files in a one archive (e.g., a zip file)? It would help me to reproduce your problem without going to every single webpage and downloading data one by one...

Sure, I probably should have done that in the first place. Here's a link to a dropbox zip folder with all the required files: https://www.dropbox.com/s/zi3nr4gs2q12hbl/maps.zip?dl=0 thanks for taking a look at this!

Please read some of my comments in the script:

library(sf)
library(raster)
library(dplyr)
library(ggplot2)

raster_data = raster("data/relief/100m/srgrhii0100a.tif")
vector_data = read_sf("data/shapefiles/coastline/Coastline.shp")

# there are two problems: -------------------------------------------------
# 1. both datasets should have the same projection, but they do not
raster_projection = projection(raster_data)
vector_projection = st_crs(vector_data)$proj4string
raster_projection
vector_projection
# 2. there are some errors in the raster object coordinates
# you can even see this error on a map at https://earthworks.stanford.edu/catalog/stanford-qh711pf3383
# (vector data is correct)
mapview::mapview(raster_data)
mapview::mapview(vector_data)

# solutions:
# 1. reproject one of the datasets
# using projectRaster (for raster data) or st_transform (for sf objects)
# read more at https://geocompr.robinlovelace.net/reproj-geo-data.html


## Hawaii
hawaii_relief <- raster_data %>%
        as("SpatialPixelsDataFrame") %>%
        as.data.frame() %>%
        rename(., value = "srgrhii0100a") %>%
        filter(value != 255)

# 2. fix the coordinates problem
# here it could be done by substracting a value from y coordinates
# I choose this value by trial-and-error
hawaii_relief$y = hawaii_relief$y - 1104000

hawaii_geo <- vector_data %>%
        st_transform(raster_projection)

ggplot(
        # define main data source
        data = hawaii_geo
) +
        # first: draw the relief
        geom_raster(
                data = hawaii_relief,
                inherit.aes = FALSE,
                aes(
                        x = x,
                        y = y,
                        alpha = value
                )
        ) +
        # use the "alpha hack" (as the "fill" aesthetic is already taken)
        scale_alpha(name = "",
                    range = c(0.4, 0.01),
                    guide = FALSE) + # suppress legend
        # add main fill aesthetic
        # use thin white stroke for municipality borders
        geom_sf(color = "red",
                fill = NA,
                size = 0.1
        )