Showing posts with label ArcPy. Show all posts
Showing posts with label ArcPy. Show all posts

Friday, August 19, 2016

Panda Dataframe as a Process Tracker (postgres example)

Sometimes you need to keep track of the number of rows processed for a given table.

Let's assume you are working in postgres and you want want to do row by row operations to do some sort of data manipulation.  Your user requires you to keep track of each row's changes and wants to know the number of failures with the updates and the number of successful updates. The output must be in a text file with pretty formatting.

There are many ways to accomplish this task, but let's use Pandas, arcpy.da Update Cursor, and some sql.


#--------------------------------------------------------------------------
def create_tracking_table(sde, tables):
    """
    creates a panadas dataframe from a sql statement
    Input:
       sde - sde connection file
       tables - name of the table to get the counts for
    Ouput:
       Panda Dataframe with column names: Table_Name, Total_Rows and
       Processed
    """
    desc = arcpy.Describe(sde)
    connectionProperties = desc.connectionProperties
    username = connectionProperties.user
    sql = """SELECT
       nspname AS schemaname,relname,reltuples
    FROM pg_class C
     LEFT JOIN pg_namespace N ON (N.oid = C.relnamespace)
    WHERE
       nspname NOT IN ('pg_catalog', 'information_schema') AND
       relkind='r' AND
       nspname='{schema}' AND
       relname in ({tables})
    ORDER BY reltuples DESC;""".format(
                                   schema=username,
                                   tables=",".join(["'%s'" % t for t in tables])
                               )
    columns = ['schemaname','Table_Name','Total_Rows']

    con = arcpy.ArcSDESQLExecute(sde)
    rows = con.execute(sql)
    count_df = pd.DataFrame.from_records(rows, columns=columns)
    del count_df['schemaname']
    count_df['Processed'] = 0

    count_df['Errors'] = 0
    return count_df



Now we have a function that will return a dataframe object from a SQL statement.  It contains 3 fields; Table_Name, Total_Rows, and Processed.  Table_name is the name of the table in the database.  Total_Rows is the length of the table.  Processed is where you are going to modify every a row gets updated successfully.  Errors is the numeric column where if an update fails, the value will be added to.

So let's use what we just made:

count_df = create_tracking_table(sde, tables)
for table in tables:
   with arcpy.da.UpdateCursor(table, "*") as urows:
      for urow in urows:
         try:
            urow[3] += 1
            urows.updateRow(urow)
            df.loc[df['Table_Name'] == '%s' % table, 'Processed'] += 1
         except:
            df.loc[df['Table_Name'] == '%s' % table, 'Errors'] += 1

The pseudo code above shows that whenever an exception is raised, 'Errors' get 1 added to it, and when it successfully updates a row 'Processed' gets updated.

The third part of the task was to output the count table to a text file which can be done easily using the to_string() method.

with open(, 'w') as writer:
   writer.write(count_df.to_string(index=False, col_space=12, justify='left'))
   writer.flush()

So there you have it.  We have a nice human readable output table in a text file.

Enjoy

Wednesday, August 3, 2016

More on Pandas Data Loading with ArcGIS (Another Example)

Large datasets can be a major problem with systems that are running 32-bit Python because there is an upper limit on memory use: 2 GB.  Most times programs fail before they even hit the 2 GB mark, but there it is.

When working with large data that cannot fit into the 2 GB of RAM, how can we push the data into DataFrames?

One way is to chunk it into groups:

#--------------------------------------------------------------------------
def grouper_it(n, iterable):
    """
    creates chunks of cursor row objects to make the memory
    footprint more manageable
    """
    it = iter(iterable)
    while True:
        chunk_it = itertools.islice(it, n)
        try:
            first_el = next(chunk_it)
        except StopIteration:
            return
        yield itertools.chain((first_el,), chunk_it) 


This code takes an iterable object (has next() defined at Python 2.7 or __next__() for Python 3.4) and makes other iterators of size n where n is a whole number (integer).

Example Usage:

import itertools
import os
import json
import arcpy
import pandas as pd


with arcpy.da.SearchCursor(fc, ["Field1", "Field2"]) as rows:
     groups = grouper_it(n=50000, iterable=rows)
     for group in groups:
         df = pd.DataFrame.from_records(group, columns=rows.fields)
         df['Field1'] = "Another Value"
         df.to_csv(r"\\sever\test.csv", mode='a')
         del group
         del df
     del groups

This is one way to manage your memory footprint by loading records in smaller bits.

Some considerations on 'n'.  I found the following effects the size of 'n': number of columns, field length, and data types.




Wednesday, July 27, 2016

Reading Spatial Data Into a Pandas Dataframe

At 10.4.x scipy is included in your basic python install, which is great!

Working with Pandas DataFrame can make life easy, especially if you need to do it quickly.


import arcpy
import pandas as pd
import sys
#--------------------------------------------------------------------------
def trace():
    """
        trace finds the line, the filename
        and error message and returns it
        to the user
    """
    import traceback
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    # script name + line number
    line = tbinfo.split(", ")[1]
    # Get Python syntax error
    #
    synerror = traceback.format_exc().splitlines()[-1]
    return line, __file__, synerror

with arcpy.da.SearchCursor(r"d:\temp\scratch.gdb\INCIDENTS_points",
                           ["OBJECTID", "SHAPE@X", "SHAPE@Y"]) as rows:
    try:
        df = pd.DataFrame.from_records(data=rows,
                                       index=None,
                                       exclude=None,
                                      columns=rows.fields,
                                      coerce_float=True)
        print ((df.columns[1], df.columns[2]))
        print ((df[df.columns[1]].mean(), df[df.columns[2]].mean()))


    except:
        print trace()




Like normal, you create an arcpy.da cursor, then pass that generator into the DataFrame's from_records().  Once the data is loaded, like in my example, you can perform operations on the frame itself.  For example let's say you needed the mean location of points.  This can be quickly done by loading in all the location XY columns (SHAPE@X and SHAPE@Y) and performing a mean call on each column.

With this method you can't control the chunksize when loading the data, so be careful of your memory.

Monday, August 31, 2015

The Walking Metadata (hermes example)

The hermes Python package was released (v1) last week, and now I'm going to show how to read the values instead a given metadata file from a feature class.

import hermes
if __name__ == "__main__":
    fc = r"c:\data_and_maps\usa\census\blkgrp.sdc\blkgrp"
    paperwork = hermes.Paperwork(dataset=fc)
    data = paperwork.convert()
    for k,v in data.iteritems():
        print k, v

In this example, I access the dataset, in this case a feature class.  To convert the metadata to a dictionary, call the convert().  Then you can use the dictionary iterator object to walk the data.

This example is just going to show something like:

metadata {some dictionary value}


Enjoy!

Friday, August 28, 2015

Work with metadata check out hermes

So you have to work with metadata? Or you love paperwork, or you just hate using the Python XML packages? Give hermes a try. It is named after And inspired by Hermes Conrad from Futurama. It is a simple set of tools to translate the Xml to a Python dictionary. Check it out here: http://github.com/Esri/hermes Fork it, follow it, improve it. I'll be posting a couple of usage examples over the next few weeks. I look forward to everyone's feedback.

Friday, July 24, 2015

A Big Thank You And ArcREST Documentation Update!

I want to thank everyone for attending my session.  It was great to hear what you are doing, and what you want to do with ArcREST, as well as just meeting all of you.  It was exciting to know that I was able to make your workday a little easier and hopefully you can REST a bit better knowing that this package is there to help you out!

I have a lot more work flows to hammer out and ideas to get working on now that I spoke to the users!

While flying back to my home, I pushed out new documentation for the v3.0.0 ArcREST, which can be found here (https://github.com/Esri/ArcREST/tree/master/docs)  If you clone the repository, the HTML documentation is pulled down.  So update your forks, or re-pull the code to get the latest documentation.  I must admit though, I am not the best with Sphinx and Latex, so it's bear bones.

Please feel free to contact me on the issue's page here, or submit a pull request to help make the help documentation look better.

Thank you for attending Esri's UC 2015, attending the ArcREST demo theater presentation, or just talking to me at the developer island.  It was a lot of fun to interact with everyone and see all the awesome things people are doing.

Keep on mapping!
~A


Tuesday, June 9, 2015

ArcREST Publishing Part 2 - SD Files

The easiest and quickest way to get data into AGOL/Portal is to publish from a service definition file, or SD file.  SD files are compressed files using 7-zip.  They contain all sorts of information about how your service behaves, what it's name is, the service's data, etc...  It's a complete package or all in one file for publishing services.

You can automate the creation of service using Python or AO, but I am not going to dive into that topic.  Instead I am going to show how you can easily publish the SD file using ArcREST to your site (AGS, AGOL or Portal).

This example is going to focus on solely the publishing of SD files to AGOL.

import arcrest
import json
import time
#--------------------------------------------------------------------------
if __name__ == "__main__":
    username = "someadminuser"
    password = "reallysecurepassword"
    sd_file = r"C:\temp\1\DemoPublishing.sd"
    proxy_port = None
    proxy_url = None

    #   Connect to the site
    #
    sh = arcrest.AGOLTokenSecurityHandler(username, password)
    admin = arcrest.manageorg.Administration(securityHandler=sh,
                                             initialize=True)
    content = admin.content

    uc = content.usercontent(username)
    ip = arcrest.manageorg.ItemParameter()
    ip.title = "DemoPublishing"
    ip.type = "Service Definition"
    ip.tags = "SD, publishing, example"
    ip.typeKeywords = "Data, Service, Service Definition"

    res = uc.addItem(itemParameters=ip, filePath=sd_file)
    itemId = res['id']
    # Publish the SD File
    #
    resPublish = uc.publishItem(fileType="serviceDefinition",
                                publishParameters=None, itemId=itemId)
    jobId = resPublish['services'][0]['jobId']
    serviceItemId = resPublish['services'][0]['serviceItemId']
    flUrl = resPublish['services'][0]['serviceurl']
    status = uc.status(itemId=serviceItemId, jobId=jobId)
    while status['status'].lower() == "processing":
        time.sleep(2)
        status = uc.status(itemId=serviceItemId, jobId=jobId)
    print 'service published.'
    print 'Bonus Demo - update some additional properties on a feature service'
    flAdmin = arcrest.agol.FeatureService(url=flUrl, securityHandler=sh).administration
    # Bonus 1- Enable edits Only
    updateDefinition = """{"hasStaticData":false,"capabilities":"Query,Editing,Create,Update,Delete","allowGeometryUpdates":true,"editorTrackingInfo":{"enableEditorTracking":false,"enableOwnershipAccessControl":false,"allowOthersToUpdate":true,"allowOthersToDelete":true}}"""
    print flAdmin.updateDefinition(json_dict=json.loads(updateDefinition))
    # Bonus 2- Enable edits with export data - adds 'Extract' to capabilities
    updateDefinition2 = """{"hasStaticData":false,"capabilities":"Query,Editing,Create,Update,Delete,Extract","editorTrackingInfo":{"enableEditorTracking":false,"enableOwnershipAccessControl":false,"allowOthersToUpdate":true,"allowOthersToDelete":true}}"""
    print flAdmin.updateDefinition(json_dict=json.loads(updateDefinition2))
    # Bonus 3- Enable Sync on services - Adds 'Sync' to capabilities
    updateDefinition3 = """{"hasStaticData":false,"capabilities":"Query,Editing,Create,Update,Delete,Sync,Extract","editorTrackingInfo":{"enableEditorTracking":false,"enableOwnershipAccessControl":false,"allowOthersToUpdate":true,"allowOthersToDelete":true}}"""
    print flAdmin.updateDefinition(json_dict=json.loads(updateDefinition3))
    # Bonus 4- Enable editor tracking on service - set enableEditorTracking = True
    updateDefinition4 = """{"hasStaticData":false,"capabilities":"Query,Editing,Create,Update,Delete,Sync,Extract","editorTrackingInfo":{"enableEditorTracking":true,"enableOwnershipAccessControl":false,"allowOthersToUpdate":true,"allowOthersToDelete":true}}"""
    print flAdmin.updateDefinition(json_dict=json.loads(updateDefinition4))
    print 'finished!'





I added some additional information to the example here.  The basic example ends at the print 'service published' part.  The additional code below shows how to enable the editing, export data, sync or track edits.  You can do this all in call, but I wanted to show how you enable each one at a time because you may only want to enable one function over another.  Notice how I get the feature service URL from the publish() then I take that URL and create a Feature Service object from the agol.FeatureService Class.  Since my user is an administrator and owner of the service, I can access the back end function to enable the additional properties.  On the Feature Service class, there is a property called administration which takes you to the hostingservice.AdminFeatureService class which allows users to update the service definitions on a hosted feature service.

I threw a lot out there, and I hope it helps!

Previous posts on publishing:
Part 1 - Publish CSV files can be found here.

Wednesday, April 8, 2015

ArcREST - Bigger and Better

I've been working hard on ArcREST (https://github.com/esri/arcrest) these days, and the package is catching on I think.  Over the last couple of months, ArcREST has been shown at both Esri's Federal Developer Summit and at the Esri Developer Summit in California.  I truly think the effort is beginning to pay off, and users are reaping the hard work put into this package.

Simply put, this is a toolbox to manage, manipulate, and control your online GIS presence.  Portal, ArcGIS Server, and ArcGIS Online can now be controlled through Python!  The modules present in the Python package are like hammers, nails, and screw drivers.  They provide the framework to do greater things.  By themselves, they do not seem like much, but together you can build your GIS house through the common scripting language of the science community.

I am constantly looking for feedback on this package, so sign up on github, and test it out with a work flow you need in these environments.  If you find a bug, have a question, or need help, please post it in the issue area of github.  

Thank you everyone!


Monday, December 22, 2014

Converting PDF to TIFF (10.3) using ArcPy

New at 10.3, is a handy tool many of us in the GIS world have wanted for a long time.  That is converting geo-referenced PDFs to Tiff files.
The help describes this tool as:
Exports an existing PDF file to a Tagged Image File Format (TIFF). If the PDF has georeference information, the TIFF can be a GeoTIFF. These TIFFs can be used as a source for heads-up digitizing and viewing in ArcMap. Both GeoPDF and ISO standards of georeferenced PDFs are supported. 

It is very straight forward to use:

import arcpy
arcpy.PDFToTIFF_conversion(in_pdf_file="C:/temp/sample.pdf", 
                           out_tiff_file="C:/temp/sample.tif", 
                           pdf_password="", 
                           pdf_page_number="1", 
                           pdf_map="Layers", 
                           clip_option="NO_CLIP", 
                           resolution="250", 
                           color_mode="RGB_TRUE_COLOR", 
                           tiff_compression="LZW", 
                           geotiff_tags="GEOTIFF_TAGS")



For this sample, I just created a map with one of the AGOL imagery layers in a blank ArcMap layout and exported it to PDF with the GeoReference information embedded.  This creates a TIFF that can be used for additional spatial data in the ArcMap session.

Enjoy,

A

Monday, November 17, 2014

Workflow for Adding Files in Parts for Portal or ArcGIS Online

Python 2.7.5 is not very good at handling big files.  Actually it stinks especially if you need to upload them via a multi-part post.  There are many 3rd party solutions, but I found issues with both poster and requests with the tool belt add-in.  They simply didn't work.  Luckily for us, the folks who designed the AGOL/Portal REST API provided an alternative way of uploading large files.  The method is called Add Item Part.  The documentation can be found here: http://resources.arcgis.com/en/help/arcgis-rest-api/index.html#/Add_Item_Part/02r300000094000000/.

When trying to figure out how to use this function, you also need to read the Add Item method as well to get how Add Item Part should be used.  The workflow is as follows:
  1. Call Add Item with no data (file data) and pass in multipart=True along with the filename=.extension.
  2. Break the file into parts, but they must be over 5 MBs except for the last chunk which signals an end of file. (I break everything into 50 MB chunks because I don't want to perform tons of POST calls)
  3. Call Add Item Part by using a POST.  Pass in the parameters partNum, which is a unique integer from 1 to 10,000.  You must pass that in, or each part will be overwritten.  Also, pass in the file via multi-part POST.
  4. Now call Commit function to tell the server, hey I'm done.  This is an asynchronous call, so don't forget to check the item's status before moving on, or you will not be able to update the item's properties.
    • If you look at your 'My Content' page on AGOL, you'll notice and item with everything stated as 'null' for the title, type, etc.. All you can do is delete the item right now.  This makes step 5 very important
  1. Update the item using the updateItem REST call.  Here we state the type, title, tags, etc.. all the good stuff you need to know so users can access the data.  

Great now we have the item in AGOL or Portal.  Simple right?  Not really, but it's a great way to get around Python's 2.7.5 annoying memory issue.  I haven't tested this on x64 python or python 3.4.  Hopefully some of my readers will post a comment letting me know if the memory issues with StringIO/cStringIO still exist when performing large file POSTs.


The add item part is in ArcREST today!

Get ArcREST: http://www.github.com/Esri/ArcREST
Also vote for my idea for a GUI builder for Python Add-Ins.

Monday, October 13, 2014

Searching for User Content on AGOL Using ArcREST

ArcGIS Online (AGOL) and Portal organizations are great, but managing the whole organization's content can be tough.  Here is a brief example on how you could search an site (in our case AGOL) and look at what a given user has shared publicly.

Searching is fairly straight forward, and do not require any security logins to perform.  You can see this by just going to www.arcgis.com and typing something into the search box.  For this example though, I am going to use login credentials just in case I needed to perform an administrative task down the line.
if __name__ == "__main__":
    username = "MY USERNAME"
    password = "MY PASSWORD"
    url = "http://www.arcgis.com/sharing"
    proxy_url = None#"127.0.0.1" # for fiddler
    proxy_port = None#"8888" # for fiddler
    q2 = "owner:esri"
    securityHandler = arcrest.AGOLTokenSecurityHandler(username,
                                                       password,
                                                       proxy_port=proxy_port,
                                                       proxy_url=proxy_url)    
    admin = arcrest.manageagol.Administration(url=url, 
                                              securityHandler=securityHandler,
                                              proxy_url=proxy_url,
                                              proxy_port=proxy_port)    
    res = admin.query_without_creds( q=q2, start=1, num=1) # Want to find only public item on owner:esri
    items = []
    total = int(res['total'])
    steps = int(total / 100)
    if (total % 100) > 0:
        steps += 1
    step = 1
    while step <= steps:
        res = admin.query_without_creds( q=q2, start= 1 + ((step-1)*100), num=100)
        for r in res['results']:
            print r # do something with item!
            del r
        del res
        step += 1    

The item is returned from the query_without_creds() as a dictionary.  It contains information about each item that is unique to that item.

The script isn't doing much with the items that it finds, but it does show how you can loop through the owner's items (esri).  If you could administer the owner's item, you could change those public items to private.

To get a better sense of what you cannot and can query, and how the syntax works, I would check out this (http://doc.arcgis.com/en/arcgis-online/reference/search.htm).  It shows the advanced search syntax in great detail.

You can get ArcREST here (www.github.com/esri/ArcREST) if you do not have it.

Happy searching!

Thursday, October 9, 2014

Using ArcREST to Export Hosted Feature Services

THIS POST IS FOR RELEASE ARCREST 2.0 AND WILL NOT WORK PERFECTLY FOR ARCREST RELEASE 3.0.

When you upload a feature class to ArcGIS Online (AGOL), you have the option to publish those feature classes. Now let's assume you have enabled editing on that data set. The question is how do you pull down the data so you can have an updated version in office? Well luckily for us, you have ArcREST to automate this task! ArcREST can be found here (www.github.com/Esri/ArcREST) if you do not have it.

This example will show how to automate the download of a hosted feature class as a shapefile to a centralized storage location.

The workflow is as follows:

  1. Connect to the AGOL site
  2. Export the item to a shapefile
  3. Download the item to disk
  4. Erase the exported item

import arcrest
import uuid
import os
import time
#   Inputs
#
itemID = "THE ITEM ID TO EXPORT"
username = "SOME USERNAME"
pw = "MY PASSWORD"
url = "http://www.arcgis.com/sharing"
filePath = r"c:\temp"
#   Logic
#
#  Create security handler
shAGOL = arcrest.AGOLTokenSecurityHandler(username, pw)
#  Connect to AGOL
org = arcrest.manageagol.Administration(url=url, securityHandler=shAGOL)
# Grab the user's content (items)
content = org.content  
usercontent = content.usercontent(username=username)
#  Create a export item with a random name using UUID
#  and export it to a zipfile
#
result =  usercontent.exportItem(title="%s" % uuid.uuid4().get_hex(),
                                     itemId=itemID,
                                     exportFormat=exportDataAs,
                                     exportParameters=None)

exportedItemId = result['exportItemId']
jobId = result['jobId']
exportItem = content.item(itemId=exportedItemId)
#   Ensure the item is finished exporting before downloading
#
status =  usercontent.status(itemId=exportedItemId, jobId=jobId, jobType="export")
    
while status['status'].lower() == "processing":
    time.sleep(3)
    status =  usercontent.status(itemId=exportedItemId,
                                 jobId=jobId,
                                 jobType="export")
filePath = exportItem.itemData(f="json", savePath=filePath)
#   Erase the exported item to clean up
#   AGOL
#
usercontent.deleteItems(items=exportItem.id)
del exportItem
print 'finishe!'

So what we have is a quick little script that downloads a file to the temp folder.  The key thing to look at is the while part of the code.  When perform actions like publish, generateFeatures, export, and createService, they are performed asynchronously on the system.  This means that you need to check if they are finished before performing further actions on the items.  The status method returns a dictionary of values, and you need to ensure that the value is not equal to 'processing', or any further actions like itemData() will return invalid information or data.  In the case of the itemData() you will download an empty zip file.

Enjoy and happy coding.

Monday, September 15, 2014

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

Friday, August 29, 2014

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.

Monday, August 11, 2014

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



Wednesday, June 18, 2014

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

Monday, June 16, 2014

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.


Wednesday, April 30, 2014

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

Tuesday, February 4, 2014

Add Geometry Attributes - ArcGIS 10.2.1 Tool

New at 10.2.1 is a simplified ways of adding geometry properties to a feature class.  The Add Geometry Attributes geoprocessing tool is a short cut from using the calculate field's geometry calculations.

Let's use the old way first:

    fc = r"some datasource"
    arcpy.AddField_management(fc, fieldName1, "DOUBLE", 
                              fieldPrecision, fieldScale)
    arcpy.AddField_management(fc, fieldName2, "DOUBLE", 
                              fieldPrecision, fieldScale)
 
    # Calculate centroid
    arcpy.CalculateField_management(fc, fieldName1, 
                                    "!SHAPE.CENTROID.X!",
                                    "PYTHON_9.3")
    arcpy.CalculateField_management(fc, fieldName2, 
                                    "!SHAPE.CENTROID.Y!",
                                    "PYTHON_9.3")
So here we create the field and then perform a calculation on the data set to get the X/Y centroid of the data.

At 10.2.1 you can clean this up and just do the following:

fc = "some dataset"
arcpy.AddGeometryAttributes_management(fc, "CENTROID")
With just two lines of code, we did what took multiple line previously to 10.2.1.

The Add Geometry Attributes tool is a time saver, and makes life just a tad easier.  It's a convenient tool that allows for quick addition of geometry attributes.

Check out more about this tool here.

Enjoy

Friday, December 20, 2013

Support Python GUI Support in ArcGIS

Please support my idea on having Python GUI support in ArcGIS.

http://ideas.arcgis.com/ideaView?id=087E00000004SmHIAU

I want to get this idea to over 1000 points before the beginning of 2014, and I need your help!

Thank you