[gdal-dev] define netcdf time dimension
Michael Sumner
mdsumner at gmail.com
Tue Apr 14 07:40:26 PDT 2015
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 <
Dominik.Schneider at colorado.edu> 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 <mdsumner at gmail.com>
> wrote:
>
>>
>>
>> On Tue, 14 Apr 2015 at 05:46 Dominik Schneider <
>> Dominik.Schneider at colorado.edu> 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
>>> gdal-dev at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150414/79fdbddd/attachment-0001.html>
More information about the gdal-dev
mailing list