Using ArcGIS and SQLite for ST_Geometry functions

spatial_sqlite_arcgis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At some point of time, you might like having your file geodatabase feature class or any other kind of spatial table loaded into a SQL database table. Maybe you are more comfortable using SQL than Python or your GIS GUI tools or maybe you already have written SQL snippets for processing your data using SQL ST_Geometry or SQL Server spatial functions and you just want to reuse them.\n",
    "\n",
    "In either case, even though your data don't reside within the database, you might accept the overhead in terms of data loading time and move your data into a PostGIS database or an SQL Server database to take advantage of ST_Geometry and SQL Server spatial functions, respectively. However, not all organizations would let you deploy a SQL Server database (the DBMS does cost some money) and even though PostGIS is open source, you might not be allowed to deploy it either on your machine or on an organization server."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this is the case, you might be surprised to hear that you will be able to use most of ST_Geometry functions without installing any additional software except for ArcGIS Desktop, which you will likely already have. Esri have supported SQLite databases for many years (in fact, it is present in versions 10.1+). This means you can use ArcMap to create an SQLite database, load your file geodatabase tables and feature classes and then use SQLite flavoured SQL to manipulate the database tables. As SQLite databases are file-based, no installation of DBMS software is required. Furthermore, Python 2 and 3 have built-in support for SQLite with the `sqlite3` module, so you won't need to install any additional packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you would rather like to perform these operations using Python (this is what I usually do), then there are multiple ways you could go. You could:\n",
    "\n",
    "* use `CopyFeatures` or `CopyRows` tools in `arcpy` to copy the tables into a SQLite database;\n",
    "* use `CreateFeatureClass` tool in `arcpy` to create a new SQLite database table and then insert your rows/features into the table;\n",
    "* use Python to load your file geodatabase data or any other data structure and then use SQL to dump the rows into SQLite database table.\n",
    "\n",
    "Below is some guidance on how this can be achieved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7.10 (default, May 23 2015, 09:40:32) [MSC v.1500 32 bit (Intel)]\n",
      "c:\\python27\\arcgis10.4\\python.exe\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print sys.version\n",
    "print sys.executable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a new `sqlite` database on disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C:\\GIS\\Temp\\Geodata20161214_155209.sqlite\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import arcpy\n",
    "from time import strftime\n",
    "\n",
    "sqlite_db = r'C:\\GIS\\Temp\\Geodata{}'.format(strftime(\"%Y%m%d_%H%M%S\"))\n",
    "if not os.path.exists(sqlite_db+'.sqlite'):\n",
    "    result = arcpy.CreateSQLiteDatabase_management(out_database_name=sqlite_db,spatial_type='ST_GEOMETRY')\n",
    "    db_path = result.getOutput(0) + '.sqlite'\n",
    "    print db_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading ArcGIS extension to support `ST_Geometry` functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<sqlite3.Cursor at 0x36f60a0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sqlite3\n",
    "conn = sqlite3.connect(db_path)\n",
    "conn.enable_load_extension(True)\n",
    "dll_path = r'c:\\Program Files (x86)\\ArcGIS\\Desktop10.4\\DatabaseSupport\\SQLite\\Windows32\\stgeometry_sqlite.dll'\n",
    "conn.execute('''SELECT load_extension('{dll_path}','SDE_SQL_funcs_init');'''.format(dll_path=dll_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def run(SQL):\n",
    "    return conn.execute(SQL).fetchall()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading features from a source file geodatabase feature class into a target `sqlite` database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import arcpy\n",
    "\n",
    "def load_into_sqlite(source_fc):\n",
    "    '''Load source feature class or table into a sqlite database'''\n",
    "    source_sr = arcpy.Describe(source_fc).spatialReference\n",
    "    source_basename = arcpy.Describe(source_fc).baseName\n",
    "    source_geometry = arcpy.Describe(source_fc).shapeType\n",
    "\n",
    "    target_fc = os.path.join(db_path,source_basename)\n",
    "    if arcpy.Exists(target_fc):\n",
    "        arcpy.Delete_management(target_fc)\n",
    "    arcpy.CreateFeatureclass_management(db_path,source_basename,source_geometry,\n",
    "                                        template=source_fc,spatial_reference=source_sr)\n",
    "\n",
    "    fields_to_exclude = ('shape','shape_area','shape_length','objectid')\n",
    "    fields = [f.name for f in arcpy.ListFields(source_fc) if f.name.lower() not in fields_to_exclude] + ['SHAPE@']\n",
    "    rows = [row for row in arcpy.da.SearchCursor(source_fc,fields)]\n",
    "    print rows[0][0:]\n",
    "\n",
    "    with arcpy.da.InsertCursor(target_fc,fields) as ins:\n",
    "        for row in rows:\n",
    "            ins.insertRow(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(u'05280', u'Bellingham', u'53', u'Washington', u'5305280', u'city', u'N', -99, 52179, <PointGeometry object at 0x71cb2b0[0x71cb420]>)\n",
      "(6380.614, u'Hawaii', u'15', u'Pacific', u'HI', 1108229.0, 1184688.0, 174, <Polygon object at 0x735c050[0x39b6f60]>)\n"
     ]
    }
   ],
   "source": [
    "fcs = [r'C:\\Program Files (x86)\\ArcGIS\\Desktop10.4\\TemplateData\\TemplateData.gdb\\USA\\cities',\n",
    "       r'C:\\Program Files (x86)\\ArcGIS\\Desktop10.4\\TemplateData\\TemplateData.gdb\\USA\\states']\n",
    "\n",
    "for fc in fcs:\n",
    "    load_into_sqlite(fc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing that features were loaded correctly (both attributes and geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(u'New York', u'POINT ( -73.94342984848072 40.66990319024563)'),\n",
       " (u'Los Angeles', u'POINT ( -118.4121144054204 34.11210653823429)')]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run('''SELECT CITY_NAME, ST_ASTEXT(SHAPE) FROM CITIES WHERE POP1990 > 3000000''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running various ST_Geometry functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(u'Los Angeles', 0.0),\n",
       " (u'Beverly Hills', 3.7541947312701076),\n",
       " (u'West Hollywood', 4.56409299421217),\n",
       " (u'Burbank', 11.806417329838071),\n",
       " (u'Culver City', 11.87289398647963)]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#get 5 nearest cities to Los Angeles, LA\n",
    "run('''SELECT CITY_NAME, ST_DISTANCE(SHAPE,(SELECT SHAPE FROM CITIES WHERE CITY_NAME = 'Los Angeles'),'Kilometer') AS DIST \n",
    "FROM CITIES WHERE STATE_NAME = 'California' ORDER BY DIST ASC LIMIT 5''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(u'Portland', -122.65770863346108, 45.5380903494073),\n",
       " (u'Eugene', -123.11339685912787, 44.05284814954899),\n",
       " (u'Salem', -123.0233092099748, 44.92434168712413),\n",
       " (u'Gresham', -122.43940384599637, 45.50369133416376),\n",
       " (u'Beaverton', -122.8178073561781, 45.475338021998084)]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#get 5 most populated cities in Oregon state\n",
    "run('''SELECT CITY_NAME, ST_X(SHAPE), ST_Y(SHAPE) \n",
    "FROM CITIES WHERE ST_WITHIN(SHAPE,(SELECT SHAPE FROM STATES WHERE STATE_NAME = 'Oregon')) = 1 \n",
    "ORDER BY POP1990 DESC LIMIT 5''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "conn.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An important note:\n",
    "\n",
    "I wasn't able to load ST_Geometry extension properly into a SQLite database using the default SQLite library which is stored at `\"C:\\Python27\\ArcGIS10.4\\DLLs\\sqlite3.dll\"`. As I have Anaconda installed, I have just replaced the SQLite library shipped with ArcGIS with the most recent library file stored at `\"C:\\Anaconda2\\DLLs\\sqlite3.dll\"`. More about this issue can be found at [this forum post](https://geonet.esri.com/thread/118976)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}