netcdf file mapping and dimensions issue

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

netcdf file mapping and dimensions issue

Warner, David
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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: netcdf file mapping and dimensions issue

Michael Sumner-2
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 <[hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: netcdf file mapping and dimensions issue

Warner, David
This looks promising!  Thank you!
Dave

On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <[hidden email]> 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 <[hidden email]> 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
>> [hidden email]
>> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: netcdf file mapping and dimensions issue

Warner, David
In reply to this post by Michael Sumner-2
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 <[hidden email]> 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 <[hidden email]> 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
>> [hidden email]
>> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: netcdf file mapping and dimensions issue

Michael Sumner-2
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 <[hidden email]> 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 <[hidden email]>
> 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 <[hidden email]> 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
>>> [hidden email]
>>> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: netcdf file mapping and dimensions issue

Warner, David
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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]>
>> 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 <[hidden email]> 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
>>>> [hidden email]
>>>> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...