I have been working on a global temperature anomaly map script to show
anomalies in 2x2 long/lat grid. With help from Roger Bivand and Andy South I
now have a script that makes the basic map with land surface background.
I use the fields package two.color() function to generate my color scheme.
my_scheme <- two.colors(n=256, start="blue", end="red", middle="white")
I choose this scheme so that "blue" reflects less than zero and "red"
reflects greater than 0 anomalies, with "white" for near zero values.
The overall map works as I hoped except that I can not figure out how to
handle missing data cells. I'd like them to show up as "grey". Currently
they are showing up as white along with the real near zero values.
How can I add a missing values color of grey to the map? Here's my overall
script.
# Load map related libraries
library(fields) # needed for image.plot()
library(sp) # needed for coordinates, SpatialPointsDataFrame
library(maptools) # neeed for wrld_simpl map
# Get world shape for map background
data(wrld_simpl) ## from maptools package
shp <- wrld_simpl
# Read source data file
link <- "http://processtrends.com/Files/global_2x2_lota_data_latest.csv";
rdf <- read.table(link, skip = 1, sep = ",", header=T)
names(rdf) <- c("i", "j", "lon", "lat", "anom")
###############################
# Convert all anom data with 9999.0000 to NA
rdf$anom[rdf$anom==9999.0000] <- NA # convert all 9999.0000 to NA
## Promote to SpatialPointsDataFrame
points_df <- rdf # make copy of original file
coordinates(points_df) = c("lon", "lat") # convert to sp file with
lon/lat cordinates
## Promote to SpatialPixelsDataFrame
pixel_df <- points_df
gridded(pixel_df) <- TRUE
## Promote to SpatialGridDataFrame
rdf_sp = as(pixel_df, "SpatialGridDataFrame")
## Establish color scheme
my_scheme <- two.colors(n=256, start="blue", end="red", middle="white")
## Generate Plot
par(mar=c(2.5,2,1,1)) ; par(oma=c(0,0,0,0)) # set plot par(mar=)
g_plot <- as.image.SpatialGridDataFrame(rdf_sp["anom"])
image.plot(g_plot, col=my_scheme,las=1,
axes=T, horizontal=T, ylim=c(-91,91), xlim=c(-181,181), zlim =
c(-6,6),
legend.mar=4,legend.shrink=0.5)
plot(shp, border="black", add=T)
--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/How-to-Show-Missing-Data-with-fields-package-two-colors-function-tp5510840p5510840.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
