R

Strsplit multiple spaces

strsplit(gsub('%', '', header), split=" +", perl=TRUE)

Draw simple map

library(maps)

library(mapdata)

library(maptools) # chama o sp e o foreign

#library(grid)

str=map(database="worldHires", plot=FALSE)

names(str)

regions=str$names

region='New Zealand(.*)Island'

(regstr=regions[grep(region, regions, ignore.case = FALSE, perl = TRUE, fixed=FALSE)])

regstr=c("New Zealand:North Island", "New Zealand:South Island")

map(database="worldHires", region=regstr,col=1)

map.axes()

image.plot legend types

imageplot_legenda.pdf

library(fields)

library(maps)

library(mapdata)

library(RColorBrewer)

library(classInt)

brks = c(0,40,200,1000,2000,4000,10000,40000)/1000.; length(brks) # nclr + 1

brksnames = c("Dense fog", "Thick fog", "Fog", "Mist", "Weak visibility.", "Moderate visibility.", "Good visibility.")

length(brksnames)

(pal = rev(c(brewer.pal(3, "Blues"), "green", "yellow", "orange", "red")))

(nclr = length(pal))

(class = classIntervals(l1$visib, n=nclr, style="fixed", fixedBreaks=brks))

(colcode = findColours(class, pal))

par(mar=c(5,4,4,10))

poly.image(xlon,xlat,l1$visib, zlim=range(class$brks),

col=attr(colcode, "palette"), breaks=class$brks)

map(database="worldHires", region=c("Portugal", "Spain"), add=TRUE, col="grey")

(atx = seq(brks[1], brks[length(brks)],length.out=length(brks)))

(atxmid = ( atx[2:length(atx)] + atx[1:(length(atx)-1)])/2)

image.plot(legend.only=TRUE, zlim=range(class$brks), col=attr(colcode, "palette"), axis.args=list(at=atxmid, labels=brksnames), legend.mar=10)

Add contours to image.plot with WRF 

(credit to Doug Nychka from the fields package)

library(RNetCDF)

library(fields)

library(maps)

library(mapdata)

filewrf = "wrfout_d01_2009-12-22_00:00:00"

xlab = expression(paste("Longitude (", degree, "E)", sep=""))

ylab = expression(paste("Latitude  (", degree, "N)", sep=""))

ncid = open.nc(filewrf)

Times = as.POSIXct(strptime(var.get.nc(ncid,"Times"), "%Y-%m-%d_%H:%M:%S"))

xlon = var.get.nc(ncid, "XLONG", start=c(1,1,1), count=c(NA,NA,1))

xlat = var.get.nc(ncid, "XLAT", start=c(1,1,1), count=c(NA,NA,1))

hgt  = var.get.nc(ncid, "HGT", start=c(1,1,1), count=c(NA,NA,1))

lmask = var.get.nc(ncid, "LANDMASK", start=c(1,1,1), count=c(NA,NA,1))

close.nc(ncid)

print(dim(hgt))

zlim = range(pretty(hgt))

print(zlim)

# interpolate data and then evaluate on a regular grid. 

bigX = cbind( c( xlon), c(xlat))

# fit interpolating surface theta is a radial distance that should be set so that 

# each point includes at least 5 and better to have 10 nearest neighbors. 

# the spacing is about .1 degrees in one direction 

# I took a factor of 4 of the spacing. 

theta.test = (xlon[33,35] - xlon[32,35])

theta.test

# [1] 0.1082153

out = fastTps(bigX, c(hgt), theta= 4*theta.test)

# evaluate interpolated surface on a grid of 120X120 points

# these will be regularly spaced in lon and lat  _not_ like the original grid

# the whole point is that contour only works in R on regularly spaced grids. 

# 

#out.p = predict.surface( out, nx=120, ny=120)

out.p = predict.surface( out, nx=60, ny=90)

png("hgt_d01.png")

hgt[lmask==0]=NA

image.plot(xlon,xlat,hgt, col=terrain.colors(64), zlim=zlim, xlab=xlab, ylab=ylab, legend.args=list(text="ter (m)",side=3, line=1))

contour(out.p, col="grey50", add=TRUE, levels=seq(200,1500,100))

map(database="worldHires", region=c("Portugal", "Spain"), add=TRUE)

map.scale(ratio=FALSE)

dev.off()

R colors

http://research.stowers-institute.org/efg/R/Color/Chart/

> colors()[grep("yellow",colors())]

[1] "greenyellow"          "lightgoldenrodyellow" "lightyellow"

[4] "lightyellow1"         "lightyellow2"         "lightyellow3"

[7] "lightyellow4"         "yellow"               "yellow1"

[10] "yellow2"              "yellow3"              "yellow4"

[13] "yellowgreen"

see his script: ColorChart.R

Plot point observations with continuous color gradient. 

The following script was greatly inspired by the Advanced Geographic Data Analysis Course of Univ. Oregon, ex4, simple maps section. Main modifications are: 

• • lateral legend created with library fields

  • color pallete with blues for negative values and reds for positive ones.

Here's the script:

library(maptools) #calls sp and foreign

library(RColorBrewer)

library(classInt)

library(fields)

# read and attach climate data

orstationc = read.csv("orstationc.csv")

print(dim(orstationc))

print(colnames(orstationc))

attach(orstationc)

# read a county outline shapefile

orotl.shp = readShapeLines("orotl.shp", proj4string=CRS("+proj=longlat"))

PlotVar = function(plotvar,class, plotclr, filename) {

png(paste("ex4_simplemaps", filename, "%02d.png", sep="_"), width=600, height=500)

colcode = findColours(class, plotclr)

cutpts = round(class$brks, digits=1)

## plot1

plot(orotl.shp)

points(lon,lat,pch=16, col=colcode, cex=2)

points(lon,lat, cex=2) # black line around each plot

par(oma=c(0,0,0,1))

legend("bottomright", legend=names(attr(colcode, "table")), fill=attr(colcode,"palette"), cex=0.7)

title(plottitle)

## plot2

op1= par(oma=c(0,0,0,4))

plot(orotl.shp, axes=TRUE)

points(lon,lat,pch=16, col=colcode, cex=2)

points(lon,lat, cex=2) # black line around each plot

title(plottitle)

op2= par(oma=c(0,0,0,1))

image.plot(legend.only=TRUE, zlim=range(cutpts), breaks=cutpts, col=attr(colcode,"palette"), axis.args=list(cex=0.5))

par(op1)

dev.off()

}

##### Precipitacao ####

arr.varname = c("pjan", "pjul", "pann")

arr.title = c("January Precipitation", "July Precipitation", "Annual Precipitation")

nclr = 8

for (i in 1:length(arr.varname)) {

varname = arr.varname[i]

plottitle = arr.title[i]

plotvar = orstationc[,varname]

# quantile

class_style = "quantile"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

plotclr = brewer.pal(nclr, "PuBu")

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

# equal

class_style = "equal"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

plotclr = brewer.pal(nclr, "PuBu")

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

# pretty

class_style = "pretty"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

plotclr = brewer.pal(nclr, "PuBu")

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

}

#### Temperature ####

arr.varname = c("tjan", "tjul", "tann")

arr.title = c("January Temperature", "July Temperature", "Annual Temperature")

nclr = 8

for (i in 1:length(arr.varname)) {

varname = arr.varname[i]

plottile = arr.title[i]

plotvar = orstationc[,varname]

# quantile

class_style = "quantile"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

if (sum(class$brks<0) > 0 && sum(class$brks>0) > 0) {

palpos = brewer.pal(sum(class$brks>0), "Reds")

palneg = brewer.pal(sum(class$brks<0), "Blues")

plotclr = c(rev(palneg),palpos)

} else {

if (sum(class$brks<0) == 0) {

plotclr = brewer.pal(nclr, "Reds")

} else {

plotclr = brewer.pal(nclr, "Blues")

}

}

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

# equal

class_style = "equal"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

if (sum(class$brks<0) > 0 && sum(class$brks>0) > 0) {

palpos = brewer.pal(sum(class$brks>0), "Reds")

palneg = brewer.pal(sum(class$brks<0), "Blues")

plotclr = c(rev(palneg),palpos)

} else {

if (sum(class$brks<0) == 0) {

plotclr = brewer.pal(nclr, "Reds")

} else {

plotclr = brewer.pal(nclr, "Blues")

}

}

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

# pretty

class_style = "pretty"

print(paste("plot", varname, class_style))

class = classIntervals(plotvar, nclr, style=class_style)

if (sum(class$brks<0) > 0 && sum(class$brks>0) > 0) {

palpos = brewer.pal(sum(class$brks>0), "Reds")

palneg = brewer.pal(sum(class$brks<0), "Blues")

plotclr = c(rev(palneg),palpos)

} else {

if (sum(class$brks<0) == 0) {

plotclr = brewer.pal(nclr, "Reds")

} else {

plotclr = brewer.pal(nclr, "Blues")

}

}

PlotVar(plotvar,class, plotclr, filename=paste(varname, class_style, sep="_"))

}