Using ArcPy to Determine the Underlying Database

Sometime you need to you what the underlying database is because you do not know what the database is.  I know it sounds strange, but it does happen.  I created a small script to solve this problem.

import arcpy
def checkDBType(dbPathOrConnectionFile):
    """ checks the db to get the database type """

    if dbPathOrConnectionFile.lower().endswith(".mdb"):
        return "PGDB"
    elif dbPathOrConnectionFile.lower().endswith(".gdb"):
        return "FGDB"
    elif dbPathOrConnectionFile.lower().endswith(".sde"):
        queries = {
            "Informix" : """select first 1 dbinfo("version", "full") from systables;""",
            "MSSQLServer" : """SELECT @@VERSION""",
            "Oracle" : """select * from v$version""",
            "PostGreSQL" : """select version()::varchar(255);"""
            }
        fail = True
        conn = arcpy.ArcSDESQLExecute(dbPathOrConnectionFile)
        for k,v in queries.iteritems():
            try:
                conn = arcpy.ArcSDESQLExecute(dbPathOrConnectionFile)
                conn.execute(v)
                return k
            except:
                pass
    return "Unknown"

Basically script just a tries the SQL snippet and it it works, then that is the database you are using. To run the sql statement, you use the arcpy.ArcSDESQLExecute()'s connection object which has a function called execute().

The results seem to be promising, but I didn't have an Informix DB to test this on. So if you have one, please let me know if that statement will work! Oracle, SQL Server and PostGreSQL all work well.

 Enjoy

Speed Up Polygon Grid Generation

There are many tools that generate regular uniform grids with geoprocessing, but I found a novel was of generating grids pretty quickly for large areas.  Rasters are fast and can be created on the fly easily using numpy.  To make the raster a feature class, ArcPy has a Raster to Polygon tool that should be available at all levels of desktop.

import arcpy
import numpy
import math
fc = r"c:\temp\somedata.shp" # assume projection is projected coordinate system not Lat/Long where measurements are in decimal degrees
width = 10000 # meters
height = 10000 # meters
#  Gets the feature classes' extent
#
desc = arcpy.Describe(fc)
extent = desc.extent
#   Set the Mask
#
arcpy.env.mask = fc
#   Calculate the number of columns
#

dist_x = math.sqrt(math.pow((extent.XMax - extent.XMin), 2))
dist_y = math.sqrt(math.pow((extent.YMax - extent.YMin),2))
col_x = int(dist_x / width)
col_y = int(dist_y / height)
if col_x % distanceMeters > 0:
   col_x += 1
if col_y % distanceMeters > 0:
  col_y += 1
#   Construct the structure numpy array with random values
#
numpyArray = (numpy.random.permutation(col_x * col_y) + 1).reshape(col_y, col_x)
#   Convert array to raster 
#
myData = arcpy.NumPyArrayToRaster(numpyArray,
                                  arcpy.Point(extent.XMin, extent.YMin),
                                  width, height)
#   Convert the raster to a polygon feature class
#
grid = r"%s\gridData" % arcpy.env.scratchGDB
grid = arcpy.RasterToPolygon_conversion(in_raster=myData,
                                            out_polygon_features=grid,
                                            simplify="NO_SIMPLIFY",
                                            raster_field="VALUE")[0]

The code takes the extent of the data and creates a grid. The mask environmental parameter essentially clips the raster to the shape of the feature class. This way you only gets cells within the given polygon area. If you wanted cells to fall outside the feature classes' shape, then just comment out that line.

So that is that.


If you haven't done so, check out ArcREST.
Also vote for my idea for a GUI builder for Python Add-Ins.

Checkout ArcREST!

Need to work with ArcGIS Online (AGOL) Feature Services?
Tired of republishing all your data to AGOL?
Want an easier way to work with REST endpoints of ArcGIS Server using Python?

Look no further, check out ArcREST on github!

Version 1 is pretty solid, and version 2 can be downloaded on the 'dev' branch.

Get it here: http://www.github.com/Esri/ArcREST

Finding the Location of Your Python File

Another issue I have run into depending on the OS, is the ability to find where the python file is located that is being run.  I like the use the __file__ property, but it's not always present.  Thanks to google, I found this:

__file__ is the path name of the file from which the module was loaded, if it was loaded from a file. The __file__ attribute is not present for C modules that are statically linked into the interpreter; for extension modules loaded dynamically from a shared library, it is the path name of the shared library file.

 This means you can't always rely on os.path.dirname(__file__) to find the folder of where your python file is located.  Here is a function that has four potential checks to get the python file's folder location:

import os
def getPythonFileLocation():
    """  returns the location of where the python file is located """
    if os.path.dirname(__file__) != "":
        return os.path.dirname(__file__)
    elif os.path.dirname(os.path.abspath(__file__)) != "":
        return os.path.dirname(os.path.abspath(__file__))
    elif os.path.dirname(os.getcwd()) != "":
        return os.path.dirname(os.getcwd())
    else:
        from inspect import getsourcefile
        return os.path.dirname(os.path.abspath(getsourcefile(lambda _:None)))

The function first tries the __file__ property using the os.path.dirname() then it tries another flavor of __file__ that I found from good old stack exchange. The 3rd check uses the os.getcwd() and finally it falls back on the inspect library.

Hopefully this will help prevent 'NoneType' issues when trying to join string paths.

Enjoy

Validating a File GeoDatabase Function

I recently created a bunch of data loading scripts that load data from some URL and copy it down local.  Nothing too special, but I noticed and interesting issue.  The file geodatabase (fgdb) would sometime randomly get corrupted and instead of the ArcGIS software seeing the fgdb as a database, it would see it as a folder in the ArcGIS software.  I create a little function to check the fgdb using the arcpy.Describe() to check the workspace ID and if it was a folder instead of a fgdb, I would erase and recreate the corrupted fgdb.  This seems to resolve the issue.

import os
import arcpy
from arcpy import env
def validate_fgdb(fgdb=env.scratchGDB):
    """
       Checks to see if the FGDB isn't corrupt.  If
       the FGDB is corrupt, it then erases the database
       and replaces it with a new FGDB of the same name.
       Input:
          fgdb - string - full path of the file geodatabase
       Output:
          Boolean
    """
    try:
        fgdb = str(fgdb)
        if os.path.isdir(fgdb) and \
           fgdb.lower().endswith(".gdb"):
            print 'workspace possible fgdb'
            desc = arcpy.Describe(fgdb)
            if hasattr(desc, "workspaceFactoryProgID") and \
               desc.workspaceFactoryProgID != "esriDataSourcesGDB.FileGDBWorkspaceFactory.1":
                shutil.rmtree(fgdb, ignore_errors=True)
                arcpy.CreateFileGDB_management(out_folder_path=os.path.dirname(fgdb),
                                               out_name=os.path.basename(fgdb))
                return True
            elif hasattr(desc, "workspaceFactoryProgID") and \
                 desc.workspaceFactoryProgID == "esriDataSourcesGDB.FileGDBWorkspaceFactory.1":
                return True
            else:
                return False

        elif os.path.isdir(fgdb) == False and \
             fgdb.lower().endswith(".gdb"):
            arcpy.CreateFileGDB_management(out_folder_path=os.path.dirname(fgdb),
                                           out_name=os.path.basename(fgdb))
            return True
        else:
            return False
    except:
        arcpy.AddError(arcpy.GetMessages(2))

Basically all it does is take a path to anything and sees if it 1 exists, 2, if it is a fgdb. I did add some additional benefits. If the path ends in ".fgdb" but does not exists, it will create the fgdb.

I have a bunch of theories on why the fgdb gets corrupted, but the leading offender is erasing tables from the fgdb causes something to get hosed in the database.  Table overwrites tend to cause this issue to happen as well.  .


Anyway enjoy and hopefully you can not eliminate this issue with your scheduled tasks or services.

Playing with GeoPackages - Shrinking the Database

Awhile back Esri announced the support of the GeoPackage, which is just a sqlite database that holds vector data, and soon raster data too.

When you add and remove data to the sqlite database, it grows but never shrinks.  Think of it like a river that rises and falls, but the water line always remains.  This waterline is the space used on your hard drive.

To reduce the footprint of the database, you can run the VACUUM command to shrink the data base.

import sqlite3
import arcpy
SQLITE_FILE = r"c:\temp\example.gpkg"
conn = sqlite3.connect(SQLITE_FILE)
conn.execute("VACUUM")
conn.close()

Pretty simple to reduce the size of the database.  What VACUUM does is reconstruct the database from scratch. This will leave the database with an empty free-list and a file that is minimal in size. Note, however, that the VACUUM can take some time to run and it can use up to twice as much temporary disk space as the original file while it is running.

Enjoy

10.2.x - Getting Data Extents as Feature Classes (4 Methods)

Over the past couple of years, I have posted many posts on creating polygon extents from data.  Today, I am going to provide 4 ways to get the extent of a polygon.

1. Use the 'Minimum Bounding Geometry' tool to create 'ENVELOPE' geometry types with or without the 'Group Option'.  This tool is available at the 'Advanced' level, so if you don't know how to/want to use python, this is your go to option to use.  The output is a feature class.
2. For whole feature class extent use the arcpy.Describe(), and write it out to disk
   

   import arcpy
   from arcpy import env
   env.overwriteOutput = True
   fc = r"c:\States.shp"desc = arcpy.Describe(fc)
   extent = desc.extent
   pts = [arcpy.Point(extent.XMin, extent.YMin),
          arcpy.Point(extent.XMax, extent.YMin),
          arcpy.Point(extent.XMax, extent.YMax),
          arcpy.Point(extent.XMin, extent.YMax)]
   array = arcpy.Array(items=pts)
   poly = arcpy.Polygon(array)
   
   arcpy.CopyFeatures_management(poly, r"%scratchgdb%\way1")
   

3. Method 1 for individual features: Use a cursor with extent XMin, YMin, XMax, and YMax, and arcpy.CreateFeatureClass()
   

   import arcpy
   from arcpy import env
   env.overwriteOutput = True
   fc = r"c:\States.shp"
   out_fc = r"%scratchgdb%\way2"
   if not arcpy.Exists(out_fc):
       arcpy.CreateFeatureclass_management(out_path="%scratchgdb%",
                                           out_name="way2", geometry_type="POLYGON",
                                           spatial_reference=arcpy.SpatialReference(4326))
   icur = arcpy.da.InsertCursor(out_fc, "SHAPE@")
   with arcpy.da.SearchCursor(fc, "SHAPE@") as rows:
       for row in rows:
           extent = row[0].extent
           pts = [arcpy.Point(extent.XMin, extent.YMin),
                  arcpy.Point(extent.XMax, extent.YMin),
                  arcpy.Point(extent.XMax, extent.YMax),
                  arcpy.Point(extent.XMin, extent.YMax)]
           array = arcpy.Array(items=pts)
           poly = arcpy.Polygon(array)
           icur.insertRow([poly])
           del array
           del poly
           del pts
           del extent
           del row
   

4. Method 2 for individual features: Use a cursor objects and arcpy.CopyFeatures()
   

   import arcpy
   from arcpy import env
   env.overwriteOutput = True
   fc = r"c:\States.shp"
   out_fc = r"%scratchgdb%\way3"
   with arcpy.da.SearchCursor(fc, "SHAPE@") as rows:
       polys = []
       array = arcpy.Array()    
       for row in rows:
           extent = row[0].extent
           array.add(extent.lowerLeft)
           array.add(extent.lowerRight)
           array.add(extent.upperRight)
           array.add(extent.upperLeft)
           array.add(extent.lowerLeft)
           polys.append(arcpy.Polygon(array))
           array.removeAll()
           del row
       del array
       arcpy.CopyFeatures_management(polys, out_fc)
       del polys
   

Here are 4 examples on how to create extents using python or the built in system tools in ArcToolbox.
 Enjoy