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

Tue May 23 17:39:49 CEST 2017

I tried a few pretty unpromising things and then pointed gdalwarp at it.
Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
it will auto-choose an output grid, and you can provide a custom one (i.e.
with a local projection.

The SDS (subdataset) naming thing is a bit awkward until you get used to
it.  This was just GDAL 2.1.2 on Ubuntu.

f <- "
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc"
download.file(f, basename(f), mode = "wb")

system('gdalwarp
HDF5:"A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/chl_ocx
chl_ocx .tif')
Creating output file that is 783P x 530L.
Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
://geophysical_data/chl_ocx.
0...10...20...30...40...50...60...70...80...90...100 - done.
r <- raster("chl_ocx ")
r
class       : RasterLayer
dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
resolution  : 0.01224719, 0.01224719  (x, y)
extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0

#plot(r)
#library(mapdata)
#map("worldHires", add = TRUE)

Cheers, Mike.

On Tue, 23 May 2017 at 23:58 Warner, David <dmwarner at usgs.gov> wrote:

> 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 <(734)%20214-9392>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]