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

Wed May 24 11:42:36 CEST 2017

Thanks a lot!  I need to fix the basics first though.  For example, even
though I can run the command and warp files at will, there is still an
issue I have yet to figure out.

The code below produces a nicely mapped file that looks like there is a
problem with the srcnodata or possibly dstnodata values as all of the pixel
values are less than zero for a scene that clearly has nonzero data.  I am
pretty stumped as I have looked at the documentation for the ocean color
data to get the srcnodata value of -32767s.  Also tried -32767, which is
shown as the fill value in SeaDAS band attributes table.  The link below
shows the result.  The data look quite like the original file as viewed in
SeaDAS but the values are clearly wrong.

r412<-raster('/Users/dmwarner/A2011202183000.L2_LAC_OC.x.nc',
var='geophysical_data/chl_ocx')
f <- "
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc"
download.file(f, basename(f), mode = "wb")
system('gdalwarp -srcnodata "None" -dstnodata None -srcnodata -32767s HDF5:"
A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/Rrs_412 output.tif')
test<-raster('/Users/dmwarner/output.tif')
plot(test)

https://cdn.rawgit.com/dmwarn/Tethys/34794e44/warprrs412.png

On Tue, May 23, 2017 at 6:04 PM, Michael Sumner <mdsumner at gmail.com> wrote:

> 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
>
>

-- 
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]]