Tuesday, 1 August 2017

Use GDAL/Python to add rasters (tiff) to an ESRI File Geodatabase


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



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...