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:
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
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:
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)
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:
My projection has a latitude_of_origin specified but no standard parallels as expected by the netcdf driver.
Best Answer
Here's some python code that does what you want, reading GDAL files that represent data at specific times and writing to a single NetCDF file that is CF-Compliant
GDAL and NetCDF4 Python can be a bit of a pain to build, but the good news is that they are part of most scientific python distributions (Python(x,y), Enthought Python Distribution, Anaconda, ...)
Update: I haven't done polar stereographic in CF-compliant NetCDF yet, but I should look something like this. Here I've assumed that
central_meridian
andlatitude_of_origin
in GDAL are the same asstraight_vertical_longitude_from_pole
andlatitude_of_projection_origin
in CF: