benfasoli
7/14/2017 - 10:22 PM
KMZ spatial data plotting in R
KMZ spatial data plotting in R
#' Generates a KMZ file
#'
#' \code{make_kmz} generates a KMZ representation of surface observations.
#'
#' @param data to be plotted on Z axis
#' @param time vector of times for data observations
#' @param lat latitude, in dd.dddd
#' @param lon longitude, in dd.dddd
#' @param filepath character path to output file, ending in '.kmz'
#' @param name to display for output kmz
#' @param cmin minimum data value to begin interpolating colors, NULL to calculate automatically
#' @param cmax maximum data value to begin interpolating colors, NULL to calculate automatically
#' @param hmin minimum interpolated height, in meters
#' @param hmax maximum interpolated height, in meters
#' @param fade points fade through array - first point is transparent, last is opaque
#' @param type 'point' plots a single point for each observation, 'line' draws a line between adjacent observations.
make_kmz <- function(data, lat, lon, filepath, time = NULL,
name = 'rkmz', species = '',
cmin = NULL, cmax = NULL,
hmin = 70, hmax = 1500,
fade = FALSE, type = 'point') {
nColor <- 64 # Number of colors to use
cbarX <- 0.06 # Colorbar dimensions (fractional)
cbarY <- 0.3
if (is.null(cmin)) cmin <- min(data, na.rm=TRUE)
if (is.null(cmax)) cmax <- max(data, na.rm=TRUE)
if (is.null(time)) time <- ''
# Assign colors to data points ------------------------------------------------
data_range <- max(data) - min(data)
height <- hmax * (data - min(data)) / data_range + hmin
cdat <- data
cdat[cdat < cmin] <- cmin
cdat[cdat > cmax] <- cmax
cidx <- ceiling((nColor - 1) * (cdat - cmin) / (cmax - cmin)) + 1
# Define desired color scale
colfun <- colorRampPalette(c('blue', 'cyan', 'green', 'yellow', 'orange', 'red'))
hex <- toupper(paste('FF', substr(colfun(nColor), 2, 7), sep='')) # in AARRGGBB
hex <- paste(substr(hex, 1, 2),
substr(hex, 7, 8),
substr(hex, 5, 6),
substr(hex, 3, 4), sep='') # in AABBGGRR
color <- hex[cidx]
if (fade){
alpha <- toupper(as.hexmode(round(seq(0, 200, length.out=length(color)))))
color <- paste(alpha, substr(color, 3, 10), sep='')
}
# ------------------------------------------------------------------------
# Set kml headers and footers
kml.header <- paste(
'<?xml version="1.0" encoding="utf-8"?>',
'<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2">',
' <Document>',
paste(' <name>',name,'</name>',sep=''),
' <Camera id="kml_cam">',
' <latitude>40.41</latitude>',
' <longitude>-111.81</longitude>',
' <altitude>20000</altitude>',
' <heading>340</heading>',
' <tilt>55</tilt>',
' <roll>0</roll>',
' <altitudeMode>absolute</altitudeMode>',
' </Camera>',
' <Style id="balloonFormat">',
' <BalloonStyle>',
' <text>',
' <![CDATA[',
' Time : $[time]',
paste(' ', name, ': $[value]'),
' ]]>',
' </text>',
' </BalloonStyle>',
' </Style>',
' <Folder>',
paste(' <name>',name,'</name>',sep=''),
paste(' <id>',name,'</id>',sep=''), '',
sep='\r\n')
kml.footer <- paste(
' </Folder>',
' <ScreenOverlay id="Colorbar">',
' <name>Colorbar</name>',
' <open>1</open>',
' <visibility>1</visibility>',
' <description/>',
' <drawOrder>0</drawOrder>',
' <color>FFFFFFFF</color>',
' <rotation>0</rotation>',
' <overlayXY x="0" xunits="fraction" y="0.5" yunits="fraction"/>',
' <screenXY x="0.01" xunits="fraction" y="0.5" yunits="fraction"/>',
paste(' <size x="',cbarX,'" xunits="fraction" y="',cbarY,'" yunits="fraction"/>',sep=''),
' <Icon>',
paste(' <href>colorbar.png</href>',sep=''),
' </Icon>',
' </ScreenOverlay>',
' </Document>',
'</kml>','',
sep='\r\n')
# ------------------------------------------------------------------------
# Generate points for kml file
if (type=='point') {
kml.pts <- c(
paste(
paste(' <Placemark id="', data, '">', sep=''),
paste(' <name>', time, data, '</name>', sep=' '),
' <styleUrl>#balloonFormat</styleUrl>',
' <visibility>1</visibility>',
' <description/>',
' <Style>',
' <LineStyle>',
paste(' <color>', color, '</color>', sep=''),
' <width>5</width>',
' </LineStyle>',
' <PolyStyle>',
' <color>00FFFFFF</color>',
' </PolyStyle>',
' </Style>',
paste(' <LineString id="', data, '">'),
' <extrude>1</extrude>',
' <tesselate>1</tesselate>',
' <altitudeMode>relativeToGround</altitudeMode>',
paste(' <coordinates>',lon, ',', lat, ',', height,' ', lon+0.00001, ',',
lat+0.00001, ',', height,'</coordinates>',sep=''),
' </LineString>',
' <ExtendedData>',
' <Data name="time">',
paste(' <value>', time, '</value>'),
' </Data>',
' <Data name="value">',
paste(' <value>', data, '</value>'),
' </Data>',
' </ExtendedData>',
' </Placemark>','',
sep='\r\n'))
} else if (type=='line') {
n <- length(lat)
lat1 <- lat[1:(n-1)]
lon1 <- lon[1:(n-1)]
height1 <- height[1:(n-1)]
lat2 <- lat[2:n]
lon2 <- lon[2:n]
height2 <- height[2:n]
kml.pts <- c(
paste(
paste(' <Placemark id="', data, '">', sep=''),
paste(' <name>', time, data, '</name>', sep=' '),
' <styleUrl>#balloonFormat</styleUrl>',
' <visibility>1</visibility>',
' <description/>',
' <Style>',
' <LineStyle>',
paste(' <color>', color, '</color>', sep=''),
' <width>5</width>',
' </LineStyle>',
' <PolyStyle>',
' <color>00FFFFFF</color>',
' </PolyStyle>',
' </Style>',
paste(' <LineString id="', data, '">'),
' <extrude>1</extrude>',
' <tesselate>1</tesselate>',
' <altitudeMode>relativeToGround</altitudeMode>',
paste(' <coordinates>',lon1, ',', lat1, ',', height1,' ', lon2, ',',
lat2, ',', height2,'</coordinates>',sep=''),
' </LineString>',
' <ExtendedData>',
' <Data name="time">',
paste(' <value>', time, '</value>'),
' </Data>',
' <Data name="value">',
paste(' <value>', data, '</value>'),
' </Data>',
' </ExtendedData>',
' </Placemark>','',
sep='\r\n'))
}
# Combine header, points, and footer
kml <- paste(kml.header, paste(kml.pts, collapse='\r\n'), kml.footer, collapse='\r\n')
# Create kml file in defined directory
dir.create(dirname(filepath),showWarnings=FALSE)
write(kml, file='.temp.kml')
# Create a colorbar
color.bar <- function(lut, cmin, cmax=-cmin, title='') {
scale = (length(lut)-1)/(cmax-cmin)
plot(c(0,10), round(c(cmin,cmax)), type='n', bty='n', xaxt='n', xlab='',
yaxt='n', ylab='', main=title, cex.axis=2, cex.main = 2)
axis(2, pretty(c(cmin, cmax), 7), las=1, cex.axis=2)
for (i in 1:(length(lut)-1)) {
y = (i-1)/scale + cmin
rect(0,y,10,y+1/scale, col=lut[i], border=NA)
}
}
png('colorbar.png', height = 500, width = 150)
color.bar(colfun(64), cmin = cmin, cmax = cmax, title = species)
dev.off()
# Convert to kmz
zip(filepath, c('.temp.kml', 'colorbar.png'))
file.remove(c('.temp.kml', 'colorbar.png'))
print(paste(filepath, ' created.'))
return(filepath)
}