DOI-USGS / gems-tools-pro

GeMS Tools for ArcGIS Pro
Creative Commons Zero v1.0 Universal
46 stars 15 forks source link

SCRIPT TO SHARE: Translate to Geopackage #54

Closed alwunder closed 1 year ago

alwunder commented 1 year ago

A Python script I wrote to translate a GeMS geodatabase to a GeoPackage. See notes in code for more info. It could probably use some error handling and cleanup, but it works as-is.

### TranslateToGeoPackage.py
### Andrew L. Wunderlich
### andrew.wunderlich@tn.gov
### 8/15/2022

### 1/20/2023  Updated to handle annotation by converting the
###    "AnnotationPolygonFeatureClass" features to points
### 2/27/2023  Ported to ArcGIS Pro GeMS toolbox

### Converts a geodatabase to an OGC GeoPackage. The script
### creates a prefix for input classes within feature datasets
### and appends it to the output class name in the GeoPackage
### so that they sort together in the output. If the input
### database is a GeMS-style database, the prefixes are
### defined (e.g., "GeologicMap" = "GM", "CrossSectionA =
### "CSA", etc.)

### The script can only convert feature classes, tables, and
### raster datasets. Other types of datasets including
### geometric network, mosaic, network, parcel fabric, raster
### catalog, schematic, terrain, tin, and topology will not be
### translated. Relationship classes are also incompatible 
### with this script.

### When adding to a toolbox:
### GENERAL:
### Name: TranslateToGeoPackage
### Label: Translate To GeoPackage
### Description: Translates a file geodatabase (gdb) to an OGC GeoPackage
### SOURCE:
### Script File: point to this file
### PARAMETERS:
### "Input database", as a Workspace (gdb to be converted)
### "Output folder" as a Folder (folder to write the GeoPackage)

import arcpy
import os
import sys
import GeMS_utilityFunctions as guf

### Set up workspace
### Use sys.argv[1] and sys.argv[2] when running script from toolbox
#inWorkspace = r"c:\temp\testing\input.gdb"
#outfolder = r"c:\temp\testing\geopackage"
inWorkspace = sys.argv[1]
outFolder = sys.argv[2]

### Set the workspace for reading the input features
arcpy.env.workspace = inWorkspace
guf.addMsgAndPrint('  Input workspace ' + inWorkspace, 0)

### Create an empty GeoPackage in the outFolder
arcpy.CreateSQLiteDatabase_management(outFolder + '/' + os.path.basename(inWorkspace)[:-4] + '.gpkg', 'GEOPACKAGE')
guf.addMsgAndPrint('  Output GeoPackage: ' + outFolder + '\\' + os.path.basename(inWorkspace)[:-4] + '.gpkg')
outWorkspace = outFolder + '/' + os.path.basename(inWorkspace)[:-4] + '.gpkg'

### Create an empty list of all input datasets (paths) that are successfully translated
translatedDatasets = []

### Create an empty list of all input datasets (paths) that do not exist in output
untranslatedDatasets = []

### Walk through input database, find ALL datasets, translate with new name based on source fds
### ALL datasets includes: Coverage, Feature, GeometricNetwork, Mosaic, Network, ParcelFabric,
### Raster, RasterCatalog, Schematic, Terrain, Tin, Table, and Topology
### This tool only converts Feature Classes, Rasters, Tables, and Annotation (Anno to points)
for root, dirs, datasets in arcpy.da.Walk(inWorkspace):
    guf.addMsgAndPrint(' ', 0)
    guf.addMsgAndPrint('  Items within feature datasets:', 0)
    for fds in dirs:
        guf.addMsgAndPrint('   Feature dataset: ' + fds, 0)
        for root2, dirs2, datasets2 in arcpy.da.Walk(fds):
            for ds2 in datasets2:
                desc = arcpy.Describe(os.path.join(root, root2, ds2))
                guf.addMsgAndPrint('    ' + desc.dataType + ' ' + os.path.join(ds2), 0)
                ### Try to translate the dataset to the target GeoPackage
                try:
                    ### Test fds name here. Add additional elif statements to test for custom GeMS feature dataset names
                    ### Use FeatureClassToFeatureClass_conversion and add fds prefix to the output name
                    if fds == 'GeologicMap':
                        fdsPfx = 'GM_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                    elif fds == 'CartographicElements':
                        fdsPfx = 'CE_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                    elif fds == 'CorrelationOfMapUnits':
                        fdsPfx = 'CMU_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                    elif fds == 'NaturalResourceMap':
                        fdsPfx = 'NRM_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                    ### CrossSection is a special case, split the fds name and use the suffix to create the Geopackage class prefix
                    elif fds.rsplit('CrossSection')[0] == '':
                        ### Get the CrossSection suffix
                        csSfx = fds.rsplit('CrossSection')[1]
                        fdsPfx = 'CS' + csSfx + '_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                    else:
                        ### If specific fds name not found, use the first 3 characters of the fds name
                        fdsPfx = fds[:3] + '_'
                        arcpy.FeatureClassToFeatureClass_conversion(ds2, outWorkspace, fdsPfx + ds2)
                        guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2), 0)
                        translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2))
                except:
                    ### This happens if the script encounters any of the incompatible class types (topology, annotation, etc.)
                    if desc.dataType == 'FeatureClass':
                        try:
                            ### Added code here for Feature to Point to handle Annotation FCs
                            ### The fdsPfx should already be set (happens just before the exception above occurs...)
                            arcpy.management.FeatureToPoint(ds2, os.path.join(outWorkspace, fdsPfx + ds2 + "_POINTS"), "CENTROID")
                            guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), fdsPfx + ds2 + "_POINTS"), 0)
                            translatedDatasets.append(os.path.join(outWorkspace, fdsPfx + ds2 + "_POINTS"))
                        except:
                            guf.addMsgAndPrint('    ...FAILED to translate ' + desc.featureType + desc.shapeType + desc.dataType + ' ' + ds2 + ' to ' + os.path.basename(outWorkspace), 1)
                            untranslatedDatasets.append(os.path.join(root, root2, ds2))
                    else:
                        guf.addMsgAndPrint('    ...FAILED to translate ' + desc.dataType + ' ' + ds2 + ' to ' + os.path.basename(outWorkspace), 1)
                        untranslatedDatasets.append(os.path.join(root, root2, ds2))
    ### Look in the root of the database for stand-alone classes
    guf.addMsgAndPrint(' ', 0)
    guf.addMsgAndPrint('   Stand-alone datasets (FCs, tables, rasters, etc.) in ' + inWorkspace + ':', 0)
    for ds in datasets:
        desc = arcpy.Describe(os.path.join(root, ds))
        guf.addMsgAndPrint('    ' + desc.dataType + ' ' + os.path.join(ds), 0)
        try:
            ### Try FeatureClass To FeatureClass conversion
            arcpy.FeatureClassToFeatureClass_conversion(ds, outWorkspace, ds)
            guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), ds), 0)
            translatedDatasets.append(os.path.join(outWorkspace, ds))
        except:
            ### NOTE: Previous version of tool ran Table try first. In that order, the script
            ### created a VAT for the raster with Value and Count fields AND copied the
            ### raster, however the VAT was always empty... Might need more testing...
            #guf.addMsgAndPrint('    ...FC to FC failed for ' + ds + '. Trying Table To Table...', 1)
            ### Try Add Raster to GeoPackage conversion
            try:
                #guf.addMsgAndPrint('    ...Attempting to translate ' + ds + ' to ' + os.path.join(os.path.basename(outWorkspace)) + ' as a raster...', 1)
                arcpy.AddRasterToGeoPackage_conversion(ds, outWorkspace, ds, "TILED")
                guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), ds), 0)
                translatedDatasets.append(os.path.join(outWorkspace, ds))
            except:
                #guf.addMsgAndPrint('    ...Add Raster to GeoPackage failed for ' + ds + '. Trying Table To Table...', 1)
                ### Try Table to Table conversion to the new GeoPackage
                try:
                    #guf.addMsgAndPrint('    ...Attempting to translate ' + ds + ' to ' + os.path.join(os.path.basename(outWorkspace)) + ' as a table...', 1)
                    arcpy.TableToTable_conversion(ds, outWorkspace, 'TAB_' + ds)
                    guf.addMsgAndPrint('    ...successfully translated to ' + os.path.join(os.path.basename(outWorkspace), 'TAB_' + ds), 0)
                    translatedDatasets.append(os.path.join(outWorkspace, 'TAB_' + ds))
                except:
                    guf.addMsgAndPrint('    ...FAILED to translate ' + desc.dataType + ' ' + ds2 + ' to ' + os.path.basename(outWorkspace), 1)
                    untranslatedDatasets.append(os.path.join(root, ds))
    break

### Show user list of unappended datasets
if len(untranslatedDatasets) > 0:
    guf.addMsgAndPrint(' ', 0)
    guf.addMsgAndPrint('Datasets that failed to translate from input database to GeoPackage ' + os.path.basename(outWorkspace) + ':', 1)
    for unappended in untranslatedDatasets:
        guf.addMsgAndPrint('   ' + unappended, 0)

### Final messages
guf.addMsgAndPrint(' ', 0)
guf.addMsgAndPrint(str(len(translatedDatasets)) + ' datasets from geodatabase ' + os.path.basename(inWorkspace) + ' translated successfully to GeoPackage ' + os.path.basename(outWorkspace) + '!', 0)
guf.addMsgAndPrint(' ', 0)