raster with hillshade in ggplot2 - episode 02

[lachlandeer]

Hi @jsta @lwasser @paleolimbot @statnmap,

moving our conversation about hillshading in #133 to a seperate issue to isolate it a little.

I've spent a little time working on this and I've found what I think is a ggplot solution that is acceptable. I think staying inside the ggplot paradigm is probably the right move - adding more and more libraries early on is probably confusing for learners. We could definitely add a tip that suggests other libraries that can do the same thing (and maybe better).

Here's the current version (after converting rasters to dataframes):

ggplot() +
    geom_raster(data = DSM_HARV_df , 
                aes(x = x, y = y, 
                     fill = HARV_dsmCrop,
                     alpha=0.8)
                ) + 
    geom_raster(data = DSM_hill_HARV_df, 
                aes(x = x, y = y, 
                  alpha = HARV_DSMhill)
                ) +
    scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) +
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0, 0.5), guide = "none") +
    ggtitle("DSM with Hillshade - NEON Harvard Forest Field Site")

And the output looks like:

which is essentially looks the same as the one in the current lesson produced with base plot up to some color softening.

If we like this approach, I'll complete the rest of the lessons with this style.

[jsta]

Looks good to me. I agree with limiting extra packages.

[lwasser]

well... i'm not sure this is doing the same thing which is why it looks fuzzy. the idea of the hillshade is the blacks represent shadows and make for a really crisp 3 dimensional looking image.

See:

Above i think the same color ramp is used for both images which is why the dimension of the image gets lost. this will be more visible if there is more terrain. here harvard is flat with mostly trees. i wonder if for now we just skip overlaying raster. OR have a breakout on doing it with baseplot?

OR show the full way which is to use the annotation_raster ? you'd do this for a publication ready plot.

i agree that adding lots of packages isn't a good idea. if we showed it with baseplot then we wouldn't have to add anything and it would look better. in this single case it could be good to just use base until a better option arises. i don't want to hold up work just call attention to the fact the plot above isn't really what we were going for originally!

but i REALLY REALLY appreicate you work on this @lachlandeer it's totally awesome that you're so invested in improving these lessons.!

[lachlandeer]

OK - let's be explicitly clear what the ggplot is doing:

1. The DSM is map is coded as follows:

# convert to a df for plotting in two steps,
# First, to a SpatialPointsDataFrame
DSM_HARV_pts <- rasterToPoints(DSM_HARV, spatial = TRUE)
# Then to a 'conventional' dataframe
DSM_HARV_df  <- data.frame(DSM_HARV_pts)

ggplot() +
    geom_raster(data = DSM_HARV_df , 
                aes(x = x, y = y, 
                     fill = HARV_dsmCrop,
                    alpha = 0.8
                     )
                ) + 
    scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) +
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0.15, 0.65), guide = "none") +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    ggtitle("DSM - NEON Harvard Forest Field Site") +
    coord_equal()

and produced this:

2. the hillshade by itself is coded as:

DSM_hill_HARV <-
  raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")

# convert to a df for plotting in two steps,
# First, to a SpatialPointsDataFrame
DSM_hill_HARV_pts <- rasterToPoints(DSM_hill_HARV, spatial = TRUE)
# Then to a 'conventional' dataframe
DSM_hill_HARV_df  <- data.frame(DSM_hill_HARV_pts)

ggplot() +
 geom_raster(data = DSM_hill_HARV_df , aes(x = x, y = y, alpha = HARV_DSMhill)) + 
    scale_alpha(range =  c(0.15, 0.65)) +
    ggtitle("Hillshade - DSM - NEON Harvard Forest Field Site") +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank())

yielding:

3. combining yields the code:

ggplot() +
    geom_raster(data = DSM_HARV_df , 
                aes(x = x, y = y, 
                     fill = HARV_dsmCrop,
                    alpha = 0.8
                     )
                ) + 
    geom_raster(data = DSM_hill_HARV_df,
                aes(x = x, y = y,
                  alpha = HARV_DSMhill)
                ) +
    scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) +
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0.15, 0.65), guide = "none") +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    ggtitle("DSM with Hillshade - NEON Harvard Forest Field Site") +
    coord_equal()

and the output:

So the main difference is that ggplot2 is only giving you three levels of grey as opposed the the ability to generate 100 with baseplot. That is where you lose the "the dimension" and gain "fuzziness." This looks like a problem with the with the easiest ggplot version. Further, I don't have a simple way to add 100 greys in ggplot, and adding an annotation_raster() feels like overkill. (as an edit: one could try the and put 100 breaks into the HARV_DSMhill via a dplyr mutate, but again feels like an extra step)

This is what is the in episode02 - PR.

I am totally OK abandoning working on this & sticking with baseplot - but then it might be awkward to switch later in the lesson once it has shapefiles. You could pitch it as 'different libraries for different problems' though, and hope no skepticism creeps in to the learners minds.

I think one thing to keep in mind is the trade off between showing the essence of the idea versus what we would want in a publication quality plot. It is not clear to me that either the ggplot or baseplot is something i would include in a paper in their current formats - yet both are instructive of a (different) approaches. Which one we go with is somewhat personal taste and I am happy to ignore mine totally 😄

[lwasser]

Great points lachlan!! Your point about showing the essence of the idea is an excellent one. What about this... what about a small breakout at the end of the lesson that discusses overlay raster plotting using raster annotation and maybe just an example using baseplot? Just so folks know other options are there... That way those folks looking for the higher quality plot get where to go. But we teach what u present which is a much more digestible approach. Thus we use all of ur work! But mention options! How does that sound? Leah Leah Wasser ******************** Director, Earth Analytics Education | Earth Lab @leahawasser https://earthdatascience.org

…

________________________________ From: Lachlan Deer <notifications@github.com> Sent: Monday, April 23, 2018 10:12:17 AM To: datacarpentry/r-raster-vector-geospatial Cc: Leah Wasser; Mention Subject: Re: [datacarpentry/r-raster-vector-geospatial] raster with hillshade in ggplot2 - episode 02 (#155) OK - let's be explicitly clear what the ggplot is doing: 1. The DSM is map is coded as follows: # convert to a df for plotting in two steps, # First, to a SpatialPointsDataFrame DSM_HARV_pts <- rasterToPoints(DSM_HARV, spatial = TRUE) # Then to a 'conventional' dataframe DSM_HARV_df <- data.frame(DSM_HARV_pts) ggplot() + geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop, alpha = 0.8 ) ) + scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) + guides(fill = guide_colorbar()) + scale_alpha(range = c(0.15, 0.65), guide = "none") + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) + ggtitle("DSM - NEON Harvard Forest Field Site") + coord_equal() and produced this: [image]<https://user-images.githubusercontent.com/4036802/39138432-2d880a1c-4720-11e8-8eb5-6c3b42bd8d45.png> 1. the hillshade by itself is coded as: DSM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif") # convert to a df for plotting in two steps, # First, to a SpatialPointsDataFrame DSM_hill_HARV_pts <- rasterToPoints(DSM_hill_HARV, spatial = TRUE) # Then to a 'conventional' dataframe DSM_hill_HARV_df <- data.frame(DSM_hill_HARV_pts) ggplot() + geom_raster(data = DSM_hill_HARV_df , aes(x = x, y = y, alpha = HARV_DSMhill)) + scale_alpha(range = c(0.15, 0.65)) + ggtitle("Hillshade - DSM - NEON Harvard Forest Field Site") + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) yielding: [image]<https://user-images.githubusercontent.com/4036802/39138311-f198445e-471f-11e8-804e-e9f1eb399193.png> 1. combining yields the code: ggplot() + geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop, alpha = 0.8 ) ) + geom_raster(data = DSM_hill_HARV_df, aes(x = x, y = y, alpha = HARV_DSMhill) ) + scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) + guides(fill = guide_colorbar()) + scale_alpha(range = c(0.15, 0.65), guide = "none") + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) + ggtitle("DSM with Hillshade - NEON Harvard Forest Field Site") + coord_equal() and the output: [image]<https://user-images.githubusercontent.com/4036802/39138559-7c325d7a-4720-11e8-9b12-f2cb2569f483.png> So the main difference is that ggplot2 is only giving you three levels of grey as opposed the the ability to generate 100 with baseplot. That is where you lose the "the dimension" and gain "fuzziness." This looks like a problem with the with the easiest ggplot version. Further, I don't have a simple way to add 100 greys in ggplot, and adding an annotation_raster() feels like overkill This is what is the in episode02 - PR. I am totally OK abandoning working on this & sticking with baseplot - but then it might be awkward to switch later in the lesson once it has shapefiles. You could pitch it as 'different libraries for different problems' though, and hope no skepticism creeps in to the learners minds. I think one thing to keep in mind is the trade off between showing the essence of the idea versus what we would want in a publication quality plot. It is not clear to me that either the ggplot or baseplot is something i would include in a paper in their current formats - yet both are instructive of a (different) approaches. Which one we go with is somewhat personal taste and I am happy to ignore mine totally 😄 — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#155 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AHS3qkhR0tSEZHMHlUP0nsHtydwfoUxjks5trf1hgaJpZM4TeinR>.

[statnmap]

I am just passing by, but if you do not want to play with alpha parameter because part of the information is difficult to see, you can also calculate your palette before the plot with cut and then add your colour vector out of the aes.
This may also help you using 100 levels of grey if you want.
I did not try it with your data but your code may looks something like:

library(dplyr)
# DSM
DSM_HARV_pts <- rasterToPoints(DSM_HARV, spatial = TRUE)
DSM_HARV_df  <- data.frame(DSM_HARV_pts)

# Hillshade
DSM_hill_HARV <-
  raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
DSM_hill_HARV_pts <- rasterToPoints(DSM_hill_HARV, spatial = TRUE)
# Calculate grey level for each tile
DSM_hill_HARV_df  <- data.frame(DSM_hill_HARV_pts) %>%
  mutate(grey = cut(layer, breaks = 10, labels = grey(seq(0,1, length.out = 10))))

# Plot
# Hillshade before as plain colour. DSM after with fixed alpha outside `aes`.
ggplot() +
    geom_raster(data = DSM_hill_HARV_df,
                aes(x = x, y = y),
                fill = DSM_hill_HARV_df$grey
                ) +
    geom_raster(data = DSM_HARV_df , 
                aes(x = x, y = y, 
                     fill = HARV_dsmCrop
                     ),
                alpha = 0.7
                ) + 
    scale_fill_gradientn(name = "Elevation", colors = rainbow(100)) +
    guides(fill = guide_colorbar())

[ErinBecker]

Hi all - chipping in late to the conversation here.

I'm not familiar with geospatial data, but from a pedagogical perspective, I would strongly recommend avoiding showing both base R and ggplot plotting (even if base R plotting is only shown in a breakaway / callout / optional exercise). From my experience, it's too easy to lose sight of the core pedagogical principles (cognitive overload, fluid representations, expert blind spot) and try to find a "perfect" solution rather than a solution that works in the context of the instructional goals for this lesson.

I fully agree with @lachlandeer that we don't need to have learners producing 100% publication quality graphics at the end of this lesson, if that goal conflicts with the goals of reducing cognitive overload, etc.

So I'm in favor (pedagogically) of moving forward with @lachlandeer's solution - although I repeat that I don't know geospatial data!

[lwasser]

sounds good to me for the time being @ErinBecker @lachlandeer