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()