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

Wed May 24 00:04:36 CEST 2017

See gdalUtils for a wrapper, or ncdump to easily get the names of variables
and other metadata. I have workers L2 here that might help you batch this
and I'm happy to help:

https://github.com/mdsumner/roc/blob/master/R/readL2.R

We want to get deeper into this soon so I'll be looking at it.

There's possibly a way to do all the sds in one GDAL call too.

Cheers, Mike

On Wed, 24 May 2017, 02:34 Warner, David <dmwarner at usgs.gov> wrote:

> Now if I can just get calls to gdal functions to work from within R....
>
>
> On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <mdsumner at gmail.com>
> wrote:
>
>>
>> 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
>>
>>
>
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]