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) 

1 Answer 1

4
+350

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.

8
  • 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? Dec 17, 2018 at 16:42
  • Do you mean this information: ggplot_build(pathPlot) ? You can get the x-y data using ggplot_build(pathPlot)$data[[2]] Dec 17, 2018 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) Dec 17, 2018 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? Dec 17, 2018 at 22:49
  • I didn't get what you mean by "...grab just the "image" part of the plot -- without all of the axes" Dec 17, 2018 at 22:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.