[R-sig-Geo] netcdf file mapping and dimensions issue

Tue May 23 15:58:17 CEST 2017

Greetings all

I am working nc files from the Ocean Biology Procesing Group.  The files
are MODIS aqua L2 ocean color files.

These files are not mapped or perhaps some would say orthorectified such
that when you create a rasterLayer from one of the contained variables and
plot it, the location and orientation of features is incorrect.  For
example, Georgian Bay, a large eastern subset of Lake Huron in the US and
Canada, may well show up west of the lake.

The dimensions of the data are read as column and row numbers rather than
latitude and longitude.

I am looking for a way to map these data in R other than my current
method.  The current method for a single nc file is shown below.  I run
this modified to fit into a parallelized foreach() loop that is used to
process folders of nc files.

The  basic steps that map the nc file data are 1) create a data.frame from
the nc files measured variable(s) and the latitude and longitude then
rasterize using vect2rast.SpatialPoints() and 2) resample the rasterLayer
created from the data.frame using resample() to an nc file I reprojected
using SeaDAS.  This seems overly complicated considering the fact that I
start with a raster and it is pretty slow (30-40 seconds per nc file in
parallel, nearly two minutes per file in the script example here).

I have two questions are 1) is there a way to eliminate the used of
vect2rast.SpatialPoints(), which is the slowest part of the process and 2)
are there any other ways to perhaps speed the process up if I can't
eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
does not provide the correct product.  It recreates what I started with,
unmapped data.  What am I missing, or is this just a function of the data I
am using?

Script is below.  I have included links to the two data files required in
the process.  One is the nc file,
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc
the other is the reprojected version of this file from SeaDAS.
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.reprojected.tif
To use them you will have to change the path to the file locations.

Thanks for your patience with my writing and poorly written code.  Any
ideas would be helpful.

Dave

library(raster)
library(plotKML)

rasterOptions(maxmemory = 1e+09)
start <- Sys.time()
#Create rasterLayers from the bands of nc file that are required for
mapping, with r412 being the measured variable of interest
r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_412')
r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_443')
r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_469')
r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_488')
r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_531')
r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_547')
r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_555')
r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_645')
r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_667')
r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var='geophysical_data/Rrs_678')

longitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var="navigation_data/longitude")
latitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
var="navigation_data/latitude")
r412  # shows that dimensions are row and column numbers, not lat and long

#Import single reprojected band of the nc file as geotif, with reprojection
done manually in SeaDAS software
correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
correct/Huron/A2011202183000.reprojected.tif")

#now rasterize and map the nc file after

d <- data.frame(x = values(longitude), y = values(latitude), r412 =
values(r412))
d$r412[is.na(d$r412)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
method='raster')
ras.412 <- raster(dd)
ras.412.crop<-crop(ras.412, myext)
#ras.412.crop[ras.412.crop<0]<-NA
ras.412.crop[ras.412.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r443 =
values(r443))
d$r443[is.na(d$r443)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
method='raster')
ras.443 <- raster(dd)
ras.443.crop<-crop(ras.443, myext)
ras.443.crop[ras.443.crop<0]<-NA
ras.443.crop[ras.443.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r469 =
values(r469))
d$r469[is.na(d$r469)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
method='raster')
ras.469 <- raster(dd)
ras.469.crop<-crop(ras.469, myext)
ras.469.crop[ras.469.crop<0]<-NA
ras.469.crop[ras.469.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r488 =
values(r488))
d$r488[is.na(d$r488)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
method='raster')
ras.488 <- raster(dd)
ras.488.crop<-crop(ras.488, myext)
ras.488.crop[ras.488.crop<0]<-NA
ras.488.crop[ras.488.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r531 =
values(r531))
d$r531[is.na(d$r531)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
method='raster')
ras.531 <- raster(dd)
ras.531.crop<-crop(ras.531, myext)
ras.531.crop[ras.531.crop<0]<-NA
ras.531.crop[ras.531.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r547 =
values(r547))
d$r547[is.na(d$r547)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
method='raster')
ras.547 <- raster(dd)
ras.547.crop<-crop(ras.547, myext)
ras.547.crop[ras.547.crop<0]<-NA
ras.547.crop[ras.547.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r555 =
values(r555))
d$r555[is.na(d$r555)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
method='raster')
ras.555 <- raster(dd)
ras.555.crop<-crop(ras.555, myext)
ras.555.crop[ras.555.crop<0]<-NA
ras.555.crop[ras.555.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r645 =
values(r645))
d$r645[is.na(d$r645)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
method='raster')
ras.645 <- raster(dd)
ras.645.crop<-crop(ras.645, myext)
ras.645.crop[ras.645.crop<0]<-NA
ras.645.crop[ras.645.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r667 =
values(r667))
d$r667[is.na(d$r667)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
method='raster')
ras.667 <- raster(dd)
ras.667.crop<-crop(ras.667, myext)
ras.667.crop[ras.667.crop<0]<-NA
ras.667.crop[ras.667.crop>90000]<- NA

d <- data.frame(x = values(longitude), y = values(latitude), r678 =
values(r678))
d$r678[is.na(d$r678)]<- 99999
coordinates(d)<-~x+y
dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
method='raster')
ras.678 <- raster(dd)
ras.678.crop<-crop(ras.678, myext)
ras.678.crop[ras.678.crop<0]<-NA
ras.678.crop[ras.678.crop>90000]<- NA

myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will cover
all of Lake Huron
correct.crop<-crop(correct, myext)
ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
End<-Sys.time()
difftime(End, start)

#plot original unmapped nc band Rrs_412 versus mapped version for gross
comparison
par(mfrow=c(2,1), mar=c(2,2,2,2))
my.orig.extent<-extent(100,475,200,600)
r412.crop<-crop(r412, my.orig.extent)
plot(r412.crop, main='Original unmapped')
plot(ras.412.r, main='Mapped in R')

-- 
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

	[[alternative HTML version deleted]]