4

I'm plotting some routes in R using the geom_path function. What I'm hoping to do is to turn the data I'm plotting into a GeoTiff (Which includes the GeoSpatial Components for Projection & Lat Long Corners) I can import into NASA WorldWind.

The artifacts I'm referencing are available here:

I made a very simple example to illustrate what I have, and what I'm trying to do:

library(rgdal)
library(ggplot2)
library(png)
library(raster)
library(tiff)


wrld <- readOGR("data" , "ne_110m_admin_0_countries")
base <- ggplot(wrld, aes(x = long, y = lat))


myDataFrame <- data.frame(Name=c("Object1","Object1","Object1","Object2","Object2","Object2"), lat=c(34,30,25,65,32,16), long=c(-118,-120,-114,-63,-108,-110)) 


route <- c(geom_path(aes(long, lat, group = myDataFrame$Name), colour = "#ffff00", size = 2, data =
myDataFrame, alpha = 0.75,
lineend = "round"))


earth <- readTIFF("HYP_LR_SR_W.tif")


pathPlot <- base +  annotation_raster(earth, -180, 180, -90, 90) + route
plot(pathPlot)

Which produces this plot:

enter image description here

The next step I'd like to do is to output the resulting plot as a GeoTIFF (which I can import into WorldWind).

I think know the commands to create a GeoTIFF in the format I want once I have a stacked raster, but I can't figure out how to wire this together to get from the route to a GeoTIFF that only includes the image itself and that includes GeoSpatial Components:

ggsave(plot=pathPlot, "pathPlot.tiff", device = "tiff")
stackedRaster <- stack("pathPlot.tiff")
xRange <- ggplot_build(pathPlot)$layout$panel_params[[1]][c("x.range")] 
yRange <- ggplot_build(pathPlot)$layout$panel_params[[1]][c("y.range")]     
extent(stackedRaster) <- extent(xRange$x.range[1],xRange$x.range[2], yRange$y.range[1],yRange$y.range[2]) 
projection(stackedRaster) <- CRS("+proj=longlat +datum=WGS84") 
writeRaster(stackedRaster, "myGeoTiff.tiff", options="PHOTOMETRIC=RGB", datatype="INT1U", overwrite=TRUE) 
mainstringargs
  • 13,563
  • 35
  • 109
  • 174

1 Answers1

4

I don't think that there is a straightforward way to coerce a ggplot object (i.e., ggproto) to a RasterStack object. I am not sure if the following solution would satisfy your requirement, but you may consider it as a workaround:

  1. Save the ggplot plot into a tiff image using ggsave
  2. Create a RasterStack object for the saved image using stack
  3. Save the RasterStack object as a GeoTiff image using writeRaster

The implementation of the above steps is as follows:

# This is your pathPlot
pathPlot <- base +  annotation_raster(earth, -180, 180, -90, 90) + route

# Remove the margins from the plot (i.e., keep the earth raster and the routes only)
pathPlot <- pathPlot + 
  theme(    
        axis.ticks=element_blank(), 
        axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.title.x=element_blank(), 
        axis.title.y=element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "null"),
        legend.position = 'none'
       ) +
       labs(x=NULL, y=NULL)

# Save the plot
ggsave(plot=pathPlot, "pathPlot.tiff", device = "tiff")

# Create a StackedRaster object from the saved plot
stackedRaster <- stack("pathPlot.tiff")

# Get the GeoSpatial Components
lat_long <- ggplot_build(pathPlot)$layout$panel_params[[1]][c("x.range","y.range")] 

# Supply GeoSpatial  data to the StackedRaster 
extent(stackedRaster) <- c(lat_long$x.range,lat_long$y.range)
projection(stackedRaster) <- CRS("+proj=longlat +datum=WGS84")

# Create the GeoTiff
writeRaster(stackedRaster, "myGeoTiff.tif", options="PHOTOMETRIC=RGB", datatype="INT1U")

Here is the produced GeoTiff image:

GeoTiff Hope it helps.

Taher A. Ghaleb
  • 5,120
  • 5
  • 31
  • 44
  • That does help to get the file itself output as a geotiff. The one big piece that is missing is the GeoSpatial part of the GeoTiff. Any thoughts on how I could extract that from the pathPlot itself or using the 'route' data structure to add that to the GeoTiff? – mainstringargs Dec 17 '18 at 16:42
  • Do you mean this information: `ggplot_build(pathPlot)` ? You can get the x-y data using `ggplot_build(pathPlot)$data[[2]]` – Taher A. Ghaleb Dec 17 '18 at 22:28
  • Ah. Yes. Great tip. I think I got the lat/lon ranges out using: xRange <- ggplot_build(pathPlot)$layout$panel_params[[1]][c("x.range")] yRange <- ggplot_build(pathPlot)$layout$panel_params[[1]][c("y.range")] extent(stackedRaster) <- extent(xRange$x.range[1],xRange$x.range[2], yRange$y.range[1],yRange$y.range[2]) projection(stackedRaster) <- CRS("+proj=longlat +datum=WGS84") writeRaster(stackedRaster, "myGeoTiff.tiff", options="PHOTOMETRIC=RGB", datatype="INT1U", overwrite=TRUE) – mainstringargs Dec 17 '18 at 22:45
  • Well that didn't paste very well. I'll add it to the main post -- I think the outstanding issue is to grab just the "image" part of the plot -- without all of the axes. Any idea how to do that? – mainstringargs Dec 17 '18 at 22:49
  • I didn't get what you mean by *"...grab just the "image" part of the plot -- without all of the axes"* – Taher A. Ghaleb Dec 17 '18 at 22:53
  • The GeoTiff has the Axes and Labels included from the Plot. I'd like to just have the Image part of the Plot (The earth raster and the lines) – mainstringargs Dec 17 '18 at 23:03
  • @mainstringargs I see what you mean. I have updated my answer. – Taher A. Ghaleb Dec 17 '18 at 23:38
  • Works great. Thank you! – mainstringargs Dec 18 '18 at 15:03