# Import a Color Composite into Google Earth (Pro)

knitr::opts_chunk\$set(eval = FALSE)

require(rgdal)
require(raster)

#Not required, probably just a left-over. TBDel in next version
#require(ggplot2)
#require(colorRamps)
#require(classInt)
#require(RStoolbox)
#source("BarryPalette.R")

rastdirAST2 <- "../../GEODATA/SUBSCENES/ASTER/ID0307566227/topcor"

Steps:

1. QGIS: Display multi-spectral image and Decide bands to be composited and stretching values
2. R: Make an rgb geotif file with the stretching you decided
3. gdal utilities: Make a pseudocolor tif version
4. gdal utilities: Make a VRT version
5. Edit the colortable in the VRT version to set the opacity of the appropriate index color to 0
6. Import VRT to Google Earth (tif must be present!)

## 1. QGIS: Display multi-spectral image in QGIS

In our example we displayed components PC1 to PC3 of a PCA image in QGIS as a RGB false color composite and interactively selected the stretching values for PC1 to PC3 as RGB. Unfortunately QGIS lacks an sliding bar for this purpose.

As you can see “para gustos los colores”: we actually found 3 different stretching sets. We use the last one in the example.

Emma

# R---min: -133.499     max: 69.9219
# G---min: -38.1393     max: 57.0268
# B---min: -25.2722     max: 21.7706

ALOBO:

# -73.2131, 69.9219
# -38.1393, 57.0268
# -25.2722, 21.7706

ALOBOv2:

# -73.2131, 69.3628
# -30.8712, 27.8837
# -13.0404, 11.486

## 2. Make an stretched rgb geotif file

We read the multi-spectral (PCA in our example) image, extract the 3 bands that were selected for the composite and apply the stretching. FInally we save this strecthed 3-bands image.

ASTbothPCA <- brick(file.path(rastdirAST2,"2008025_AST_07XT_BOTH_clip_p1_GCPOL2v2_PCA.tif"))
names(ASTbothPCA) <- paste0("PC",1:9)
summary(ASTbothPCA)

Make sure the projection is right:

projection(ASTbothPCA)

Subset the 3 first bands

rgb <- brick(subset(ASTbothPCA,1:3))

Stetching values were interactively set to:

stretchlim <- matrix(c(-73.2131, 69.3628,
-30.8712, 27.8837,
-13.0404, 11.486) ,ncol=2,byrow=TRUE)

Apply stretching

cellStats(rgb,range)
for(i in 1:3){
b <- (255 - 0)/(stretchlim[i,2] - stretchlim[i,1])
a <- 255 - b*stretchlim[i,2]
rgb[[i]] <- rgb[[i]]*b + a
}
rgb[rgb<=1] <-1; rgb[rgb>255] <-255
rgb <- round(rgb)
summary(rgb)

Write 3-bands tif image

writeRaster(rgb,file="DecepPCA1_3v2", NAflag=0,
format="GTiff",type="INT1U",overwrite=TRUE)

Note that the image file must be Unsigned integer of 1 byte (required for next step)

## 3. Convert 3-layer tif to pseudocolor tif

We use gdalutilities (http://www.gdal.org/gdal_utilities.html) that should have been installed along the rest of GDAL (http://www.gdal.org/index.html)

system("rgb2pct.py DecepPCA1_3v2.tif DecepPCApseudocolorv2.tif")

Check the output file with QGIS. Note you must dispaly without additional stretching (i.e., take min/max values for each band)

We should be ready for importing into GE here, but as, unfortunately, null values are not saved as transparent in the colortable of the pseudocolor tif: we have to find a way to edit that colortable.
Despite some references (see rgb2GE_tests_log.R), there is no way to set a given color to transparent in the pseudocolor tif itself. The workaround solution is making a VRT version and edit the colortable.

## 4. Make a VRT version

We use gdalutilities again:

system("gdal_translate -of VRT DecepPCApseudocolorv2.tif DecepPCApseudocolorv2.vrt")

Now display the vrt image (e.g. in QGIS) and check the value of the part that should become transparent. In this example we find that the value is 252.
The vrt file is just a xml file that can be opened with a simple text editor. See http://www.gdal.org/gdal_vrttut.html for details on the vrt format but all you need to know right now is that the vrt file includes the relative path and file name of the image it comes from (DecepPCApseudocolorv2.tif in this example) and the color table in text format with 256 lines. As far as I know, while the vrt image can be displayed in QGIS, colors are all set to opaque and there is no way to modify transparency for a given color in the color table. Thus we have to change opacity by editing the vrt file with a text editor.

As the color we want to make transparent is black, find the line:

# <Entry c1="0" c2="0" c3="0" c4="255" />  

and edit the c4 field to become 0:

# <Entry c1="0" c2="0" c3="0" c4="0" />

c4 codes the opacity, which we have changed from its maximum (255) to its minimum (0), thus making that given color (black) become transparent. You can display the vrt file in QGIS again to check that those parts that were black are transparent now.

## 5. Import into Google Earth Pro

GE Pro version is required, but Pro version is free now.
Note you need both the VRT and the TIF files (here DecepPCApseudocolor.vrt and DecepPCApseudocolor.tif), with the tif file in the same relative directory as it was when the vrt file was created. Best having both in the same directory.

Use File/Import: as the vrt file is fully georeferenced and with an standard projection, Google Earth will be able to put it at the correct site (provided the tif file is present in the correct directory).
I can imagine that in case of large files, it would be better providing the raster in GE’s PseudoMercator projection (EPSG: 3857). In that case, you should reproject the initial rgb brick file.