Moving from gdal-dev mailing list:
On Mon, Sep 2, 2013 at 7:09 PM, David Shean wrote:
Hi list,
I'm trying to package a timeseries of GTiff rasters with identical
projection/extent/resolution as a single NetCDF file for distribution. I've
spent the past hour consulting the online doc and playing with
gdal_translate, gdalbuildvrt and gdalwarp without any success.
Is there an easy way to do this using existing gdal command line
utilities? I figured I'd ask before resorting to a custom solution using
the NetCDF Python API.
Thanks.
-David
On Tue, Sep 3, 2013 at 10:15 AM, Etienne Tourigny wrote:
what you want is probably outside the scope of gdal. It would require some
clever metadata management so that gdal_translate puts them in a single
file...
I would advise you convert them all to netcdf using gdal_translate and then
use python-netcdf4 (not the one from numpy/scipy) to stack them in the
temporal dimension.
On Tue, Sep 3, 2013, at 7:55 AM, "Signell, Richard" wrote:
David,
If you post your question on the GIS stackexchange group
https://gis.stackexchange.com/
I will provide an example code that should be helpful.
-Rich
====================
Update 9/3/13 17:04 PDT
Here is gdalinfo output for one of my input datasets:
<pre><code>
gdalinfo 20120901T2024_align_x+22.19_y+3.68_z+14.97_warp.tif
Driver: GTiff/GeoTIFF
Files: 20120901T2024_align_x+22.19_y+3.68_z+14.97_warp.tif
Size is 10666, 13387
Coordinate System is:
PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",70],
PARAMETER["central_meridian",-45],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-211346.063781524338992,-2245136.291794800199568)
Pixel Size = (5.000000000000000,-5.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -211346.064,-2245136.292) ( 50d22'39.70"W, 69d23'55.59"N)
Lower Left ( -211346.064,-2312071.292) ( 50d13'22.38"W, 68d48'10.75"N)
Upper Right ( -158016.064,-2245136.292) ( 49d 1'33.33"W, 69d26'16.42"N)
Lower Right ( -158016.064,-2312071.292) ( 48d54'35.06"W, 68d50'27.28"N)
Center ( -184681.064,-2278603.792) ( 49d38' 1.32"W, 69d 7'17.04"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
NoData Value=-32767
</code></pre>
Following up on Luke's suggested approach.
The vrt generation works fine:
gdalbuildvrt -separate newtest.vrt *warp.tif
<VRTDataset rasterXSize="10666" rasterYSize="13387">
<SRS>PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]</SRS>
<GeoTransform> -2.1134606378152434e+05, 5.0000000000000000e+00, 0.0000000000000000e+00, -2.2451362917948002e+06, 0.0000000000000000e+00, -5.0000000000000000e+00</GeoTransform>
<VRTRasterBand dataType="Float32" band="1">
<NoDataValue>-3.27670000000000E+04</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="1">20110619T2024_align_x+15.51_y+1.15_z+12.10_warp.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="10666" RasterYSize="13387" DataType="Float32" BlockXSize="256" BlockYSize="256" />
<SrcRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<DstRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<NODATA>-32767</NODATA>
</ComplexSource>
</VRTRasterBand>
<VRTRasterBand dataType="Float32" band="2">
<NoDataValue>-3.27670000000000E+04</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="1">20110802T2024_align_x+16.33_y+2.14_z+12.02_warp.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="10666" RasterYSize="13387" DataType="Float32" BlockXSize="256" BlockYSize="256" />
<SrcRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<DstRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<NODATA>-32767</NODATA>
</ComplexSource>
</VRTRasterBand>
...
But when I attempt to translate to nc, I get the following error:
<pre><code>
gdal_translate -of netcdf newtest.vrt newtest.nc
Input file size is 10666, 13387
Warning 1: Variable has 0 dimension(s) - not supported.
0...10...20...30...40...50ERROR 1: netcdf error #-62 : NetCDF: One or more variable sizes violate format constraints .
at (netcdfdataset.cpp,SetDefineMode,1574)
ERROR 1: netcdf error #-39 : NetCDF: Operation not allowed in define mode .
at (netcdfdataset.cpp,IWriteBlock,1435)
ERROR 1: netCDF scanline write failed: NetCDF: Operation not allowed in define mode
ERROR 1: An error occured while writing a dirty block
...ERROR 1: netcdf error #-39 : NetCDF: Operation not allowed in define mode .
at (netcdfdataset.cpp,IWriteBlock,1435)
ERROR 1: netCDF scanline write failed: NetCDF: Operation not allowed in define mode
ERROR 1: netcdf error #-62 : NetCDF: One or more variable sizes violate format constraints .
at (netcdfdataset.cpp,~netCDFDataset,1548)
</code></pre>
So upon closer inspection, it appears that gdal is unhappy with the polar stereographic projection I'm using (EPSG:3413). See lines 1570-1582 of netcdfdataset.cpp:
https://code.vpac.org/gitorious/gdal-netcdf-testing/gdal-netcdf-driver/blobs/8fa3582669969ad4d55e461f5846b3ed33727f63/gdal/frmts/netcdf/netcdfdataset.cpp
My projection has a latitude_of_origin specified but no standard parallels as expected by the netcdf driver.