[gdal-dev] define netcdf time dimension
Michael Sumner
mdsumner at gmail.com
Tue Apr 14 07:41:49 PDT 2015
Sorry, edit - the setZ line should be:
r <- setZ(r, dt0)
On Wed, 15 Apr 2015 at 00:40 Michael Sumner <mdsumner at gmail.com> 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 <
> 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/b03e98da/attachment-0001.html>
More information about the gdal-dev
mailing list