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

Tue May 23 18:33:40 CEST 2017

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

	[[alternative HTML version deleted]]