0
votes

I downloaded a netcdf file that shows chlorophyll concentrations. I am trying to plot this information as a raster file over a shapefile map. When I plot the raster layer, I am having trouble with two things

  1. I can't figure out how to display coordinates with the chlorophyll layer - it won't let me turn my lat and lon into a raster layer

  2. I want to overlay the raster layer over a shapefile map ... but it seems like the NA values are also showing a a color, so it won't overlay over anything?

Here is the code I have so far:

library(sp)
library(maptools)
library(RNetCDF)
library(ncdf4)
library(raster)

# retrieving chlorophyll info
attributes(mycdf$var)
attributes(mycdf$var$chlor_a)

# get lat lon info
attributes(mycdf$dim)
nc_lat <- ncvar_get(mycdf, attributes(mycdf$dim)$names[1])
nc_lon <- ncvar_get(mycdf, attributes(mycdf$dim)$names[2])
mycdf <- nc_open(file.choose(), verbose = TRUE, write = FALSE)

rasternew <- raster(chlor)
projection(rasternew)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 
+no_defs")
rasternew<- t(rasternew)

# get map
library(rworldmap)
library(rworldxtra)
newmap<- getMap(resolution = "high")
View(newmap)
projection(newmap)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 
+no_defs")

# plot
brks <- c(0, 01, 0.2, 0.3, 0.4, 0.5)
col<- c("aliceblue", "lightblue", "green", "yellow", "orange", "red")
par(mfrow=c(1,1), font.axis = 4, mai = rep(.1,4), pin = c(5, 4))
plot(rasternew, breaks = brks, col=col, legend = F)
par(mfrow=c(1,1), new = TRUE, font.axis = 4, mai = rep(.1,4), pin = 
c(5, 4))
plot(newmap, xlim=c(150,152), ylim=c(-12, -9), asp = 1, axes = F, col = 
"lightgrey", yaxt = "n", xaxt = "n", add = T)
legend(x='bottomleft', legend = c("0", "0.1", "0.2", "0.3", "0.4", 
"0.5"), fill = col)

enter image description here

1
Thank you so much @RobertH that worked. I was just making it too complicated - Shannon Murphy

1 Answers

1
votes

Making a Raster* object from data in a ncdf file can normally be done like this:

library(raster)
r <- raster("filename.nc", var="chlor_a")
plot(r)

Then you should be able to add SpatialPolygons like this:

plot(pols, add=TRUE)