oc_sst_chl.R

u1 <- "https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/T20171522017181.L3m_MO_NSST_sst_9km.nc"
f1 <- basename(u1)

u2 <- "https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/T20171522017181.L3m_MO_CHL_chlor_a_9km.nc"
f2 <- basename(u2)

curl::curl_download(u1, f1)
curl::curl_download(u2, f2)

library(raster)

## Loading required package: sp

## note that scaling for sst happens automatically
## ncdf4 does it by default, and raster gets it
sst <- raster(f1, varname = "sst")

## Loading required namespace: ncdf4

## chl is not scaled internally so nothing changes on read
chl <- raster(f2, varname  = "chlor_a")

## crop (this breaks delay and pulls data from file)
ex <- extent(-100, -68, 12, 44)
sst <- crop(sst, ex, snap = "out")
chl <- crop(chl, ex, snap = "out")

## get some idea of the masking
plot(brick(sst >= 15 & sst <= 25,
           chl >= 2 & chl <= 5), col = c("grey", "firebrick"))

result <- (sst >= 15 & sst <= 25) &
          (chl >= 2 & chl <= 5)

## get as locations
xy <- xyFromCell(result, which(values(result)))
plot(xy, asp = 1/cos(mean(xy[,2]) * pi/180), pch = "x", cex = 0.3, col = "firebrick")
library(mapdata)

## Loading required package: maps

maps::map("worldHires", add = TRUE) ## needs mapdata pkg