I'm using GDAL python to create a bunch of tiff files from a netcdf file.
def modelout_to_geotif(netcdfile, outdir, VarName):
""" Given the path to the netcdf output,
produces geotiffs for each band of the given variable name.
Input::
netcdfile: string path to netcdf file
VarName: string variable name in netcdf file (e.g. 'GridPrec')
Output::
creates geotiffs for each 'band' (e.g. time step). Output is written
to current working directory.
"""
try:
nci = gdal.Open('NETCDF:{0}:{1}'.format(netcdfile, VarName))
ncd = Dataset(netcdfile)
except:
print "Could not open input file: {0}".format(netcdfile)
return
try:
geotransform = nci.GetGeoTransform()
except:
print "Could not get geotransform.\n"
return
try:
projection = nci.GetProjection()
except:
print "Could not get projection.\n"
return
numbands = nci.RasterCount
times = [o2d(t) for t in ncd.variables['time']]
print "Found {0} bands to process for {1}".format(numbands, VarName)
for i in range(numbands):
time = times[i]
band = nci.GetRasterBand(i+1)
raster = band.ReadAsArray()
y,x = raster.shape
output_file = "{0}_{1}.tif".format(VarName, time.strftime('%Y%m%d_%H'))
subdir = os.path.join(outdir, VarName)
if not os.path.exists(subdir):
os.mkdir(subdir)
output_file = os.path.join(subdir, output_file)
if os.path.exists(output_file):
continue
# Create gtif
driver = gdal.GetDriverByName("GTiff")
#TODO: Make this more generic
dst_ds = driver.Create(output_file, x, y, 1, gdal.GDT_CFloat32 )
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( geotransform )
# set the reference info
dst_ds.SetProjection( projection )
# write the band
dst_ds.GetRasterBand(1).WriteArray(raster)
Now I would like to add those to a mosaic catalog in an ESRI File Geodatabase. I've managed to get it working using arcpy, but it's really slow, and I'd prefer to use GDAL. Is there a way? Below is my arcpy script:
import os
try:
import arcpy
from arcpy import env
except:
raise ImportError("Unable to import arcpy - are you using ESRI Python?")
env = 'P:/Enki'
def add_to_mosaic(filegdb, raster_dir):
print env
#Add Raster Dataset type Raster to FGDB Mosaic Dataset
#Calculate Cell Size Ranges and Build Boundary
#Build Overviews for Mosaic Dataset upon the 3rd level Raster Dataset pyramid
#Apply TIFF file filter
#Build Pyramids for the source datasets
varnames = ['GridPrec']#,
#'GridSnowOut',
#'GridRad',
#'GridTemp',
#'GridSoilOut',
#'GridRunoff']
for VarName in varnames:
gdbname = filegdb
mdname = os.path.join(gdbname, VarName)
messages = 'Failed'
sr = arcpy.SpatialReference()
# sr.name = 'WGS_1984_UTM_Zone_32N'
sr.factoryCode = 32632
sr.create()
print sr
try:
noband = 1
res = None
#TODO: get from input data
print 'trying to create %s' % VarName
res = arcpy.CreateMosaicDataset_management(gdbname, mdname, sr, noband)
except:
print res.getMessages()
rastype = "Raster Dataset"
inpath = os.path.join(raster_dir, VarName)
updatecs = "UPDATE_CELL_SIZES"
updatebnd = "UPDATE_BOUNDARY"
updateovr = "UPDATE_OVERVIEWS"
maxlevel = "2"
maxcs = "#"
maxdim = "#"
spatialref = "#"
inputdatafilter = "*200701*.tif"
subfolder = "NO_SUBFOLDERS"
duplicate = "EXCLUDE_DUPLICATES"
buildpy = "BUILD_PYRAMIDS"
calcstats = "CALCULATE_STATISTICS"
buildthumb = "NO_THUMBNAILS"
comments = "Add Raster Datasets"
forcesr = "#"
res = arcpy.AddRastersToMosaicDataset_management(mdname, rastype, inpath, updatecs, \
updatebnd, updateovr, maxlevel, maxcs, \
maxdim, spatialref, inputdatafilter, \
subfolder, duplicate, buildpy, calcstats, \
buildthumb, comments, forcesr)
print res.getMessages()