Sorry, edit - the setZ line should be: r <- setZ(r, dt0)
On Wed, 15 Apr 2015 at 00:40 Michael Sumner <[email protected]> wrote: > Well, FWIW - > > Here the CRS doesn't round-trip properly, partly as the PROJ.4 string in R > is not sufficient, and the date type gets recorded as raw integer (that's > ok depending what your target environment will need). Using gdal_translate > you could augment the netcdf created here with a proper WKT version, so > maybe this is useful in a limited context. > > I tried it so you might as well see the result: > > library(raster) > r <- stack("swe_sub.mdl") > r > #class : RasterStack > #dimensions : 190, 352, 66880, 1 (nrow, ncol, ncell, nlayers) > #resolution : 30, 30 (x, y) > #extent : 444345.1, 454905.1, 4430445, 4436145 (xmin, xmax, ymin, > ymax) > #coord. ref. : +proj=utm +zone=13 +datum=NAD83 +units=m +no_defs > #names : Band.1 > > dt0 <- seq(as.POSIXct("2014-03-01 00:00:00", tz = "UTC"), by = "1 hour", > length = nlayers(r)) > r <- setZ(r, as.character(dt0)) > writeRaster(r, "swe_sub.nc") > > system("gdalinfo swe_sub.nc") > > Warning 1: dimension #2 (easting) is not a Longitude/X dimension. > Warning 1: dimension #1 (northing) is not a Latitude/Y dimension. > Driver: netCDF/Network Common Data Format > Files: swe_sub.nc > Size is 352, 190 > Coordinate System is `' > Origin = (444345.080000000016298,4436145.280000000260770) > Pixel Size = (30.000000000000000,-30.000000000000000) > Metadata: > easting#long_name=easting > easting#units=meter > layer#_FillValue=-3.4e+38 > layer#long_name=layer > layer#max=3.695244550704956 > layer#min=0 > layer#missing_value=-3.4e+38 > layer#projection=+proj=utm +zone=13 +datum=NAD83 +units=m +no_defs > layer#projection_format=PROJ.4 > NC_GLOBAL#Conventions=CF-1.4 > NC_GLOBAL#created_by=R, packages ncdf and raster (version 2.3-33) > NC_GLOBAL#date=2015-04-14 04:29:01 > NETCDF_DIM_EXTRA={time} > NETCDF_DIM_time_DEF={1,6} > NETCDF_DIM_time_VALUES=1393632000 > northing#long_name=northing > northing#units=meter > time#long_name=time > time#units=seconds since 1970-01-01 00:00:00 > Corner Coordinates: > Upper Left ( 444345.080, 4436145.280) > Lower Left ( 444345.080, 4430445.280) > Upper Right ( 454905.080, 4436145.280) > Lower Right ( 454905.080, 4430445.280) > Center ( 449625.080, 4433295.280) > Band 1 Block=352x1 Type=Float32, ColorInterp=Undefined > NoData Value=-3.39999999999999996e+38 > Metadata: > _FillValue=-3.4e+38 > long_name=layer > max=3.695244550704956 > min=0 > missing_value=-3.4e+38 > NETCDF_DIM_time=1393632000 > NETCDF_VARNAME=layer > projection=+proj=utm +zone=13 +datum=NAD83 +units=m +no_defs > projection_format=PROJ.4 > > ## context of the R session > sessionInfo() > R version 3.1.3 (2015-03-09) > Platform: x86_64-pc-linux-gnu (64-bit) > Running under: Ubuntu 14.04.2 LTS > > locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] raster_2.3-33 sp_1.0-17 > > loaded via a namespace (and not attached): > [1] grid_3.1.3 lattice_0.20-31 ncdf4_1.13 rgdal_0.9-2 > > > On Wed, 15 Apr 2015 at 00:20 Dominik Schneider < > [email protected]> wrote: > >> Hi Mike - Thanks for the response and sorry this wasn’t quite clear. IDL >> data cubes are of the ENVI format and thus also have an associated ascii >> .hdr file. A dataset with 10 bands is here: >> https://dl.dropboxusercontent.com/u/5962491/IDLENVIswe.zip >> There is no time data in the .mdl file, but I know that the 4416 bands >> are hourly timesteps from 0:00 March 1 to 23:00 Aug 31 >> >> I could do this in R but was trying to avoid it because R is always so >> frustratingly slow for these things…the gdal commandline utilities are much >> more awesome! >> >> iMac:~/Downloads $ gdalinfo swe.mdl >> Driver: ENVI/ENVI .hdr Labelled >> Files: swe.mdl >> swe.hdr >> Size is 352, 190 >> Coordinate System is: >> PROJCS["UTM_Zone_13N", >> GEOGCS["GCS_North_American_1983", >> DATUM["North_American_Datum_1983", >> SPHEROID["GRS_1980",6378137.0,298.257222101]], >> PRIMEM["Greenwich",0.0], >> UNIT["Degree",0.0174532925199433]], >> PROJECTION["Transverse_Mercator"], >> PARAMETER["False_Easting",500000.0], >> PARAMETER["False_Northing",0.0], >> PARAMETER["Central_Meridian",-105.0], >> PARAMETER["Scale_Factor",0.9996], >> PARAMETER["Latitude_Of_Origin",0.0], >> UNIT["Meter",1]] >> Origin = (444345.080000000016298,4436145.280000000260770) >> Pixel Size = (30.000000000000000,-30.000000000000000) >> Image Structure Metadata: >> INTERLEAVE=BAND >> Corner Coordinates: >> Upper Left ( 444345.080, 4436145.280) (105d39' 9.74"W, 40d 4'25.45"N) >> Lower Left ( 444345.080, 4430445.280) (105d39' 7.97"W, 40d 1'20.58"N) >> Upper Right ( 454905.080, 4436145.280) (105d31'43.92"W, 40d 4'27.72"N) >> Lower Right ( 454905.080, 4430445.280) (105d31'42.49"W, 40d 1'22.85"N) >> Center ( 449625.080, 4433295.280) (105d35'26.03"W, 40d 2'54.21"N) >> Band 1 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 2 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 3 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 4 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 5 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 6 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 7 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 8 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 9 Block=352x1 Type=Float32, ColorInterp=Undefined >> Band 10 Block=352x1 Type=Float32, ColorInterp=Undefined >> ... >> Band 4416 Block=352x1 Type=Float32, ColorInterp=Undefined >> >> >> >> Dominik Schneider >> o 303.735.6296 | c 518.956.3978 >> >> >> On Tue, Apr 14, 2015 at 12:07 AM, Michael Sumner <[email protected]> >> wrote: >> >>> >>> >>> On Tue, 14 Apr 2015 at 05:46 Dominik Schneider < >>> [email protected]> wrote: >>> >>>> Hi - >>>> >>>> I have some .mdl files from IDL with .hdr header files that I’d like to >>>> convert to netcdf. >>>> The following code produces a netcdf file with a subdataset for each >>>> band in the mdl file. Is there any way to define the bands as the time >>>> dimension, in this case 4416 hourly timesteps? >>>> >>>> Or is there another tool that can convert this? >>>> >>>> gdal_translate -of NetCDF swe.mdl swe.nc >>>> >>>> >>> Can you list the metadata printed out by gdalinfo swe.mdl? What driver >>> is used? >>> >>> Does that source file have time-metadata inside it each subdataset or do >>> you need to assign it? >>> >>> Generally, the subdatasets are not considered as a time series - they >>> are more like multiple variables for the same observation (whereas a 3rd >>> and higher dimension/s are often unrolled into a multi-attributed 2D layer >>> and the time/z step is stored on each band - it's a bit of a mix/clash of >>> worlds). >>> >>> Can you point to and example file somewhere? I imagine you'll need a >>> programmed solution, but it should be pretty easy with R or python or >>> similar. If R is of interest you can try this and take your question to >>> the R community: >>> >>> library(raster) >>> r <- stack("swe.mdl") >>> r <- setZ(r, [whatever the vector of dates should be]) ## pseudocode >>> writeRaster(r, "swe.nc") >>> >>> That all assumes quite a lot, including available NetCDF versions and >>> packages, required structure in the output, interpretation of the data read >>> in, - so it's just included as an indication how it might be done. >>> >>> Cheers, Mike. >>> >>> >>>> Thanks >>>> Dominik >>>> _______________________________________________ >>>> gdal-dev mailing list >>>> [email protected] >>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>> >>> >>
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
