2

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()
4
  • 1
    As far as I know, there is no raster support for the filegdb api released by esri. I wouldn't count on raster support any time, but who knows?
    – user10353
    Commented May 29, 2013 at 5:40
  • So you mean that the api released and accessible to gdal won't likely include the 'CreateMosaicDataset_management' functionality?
    – John
    Commented May 29, 2013 at 7:12
  • 2
    I think Kyle is right, the only way to manage Mosaic Datasets currently is through ArcPy. I've also run into performance bottleneck issues with this; it's quite frustrating since GDAL is so fast for reprojection etc.
    – greenlaw
    Commented May 29, 2013 at 14:19
  • 2
    I had started an ideas vote up thing for raster support, it garnered little support. My attempts to discuss it on arcmap forums was equally ignored. The FileGDB currently doesn't support rasters.
    – user10353
    Commented May 31, 2013 at 14:46

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.