Creating Polygons from Raster Extents

There are many like it, but this one is mine. Catch that reference?

 OK, so here is a script that will allow the user to select a folder directory. The script will scan the directory for any raster files. It will grab the extents of the raster and build a polygon feature class. When building a mosaic I like to also know the compression type, number of bands, no data value, pixel type, etc. Those options will also be added to the polygon per raster tile.

 Quick note, I wrote this with the idea that a person could add this script to their toolbox to have a nice GUI interface. The only parameter required is the input file path.



 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
# Name: Raster Extents to Polygons
# Purpose: This script will create a polygon for each raster file extent in a directory.
# Author:   Nicole Ceranek 
# ArcGIS Version    10.1
# Python Version:   2.7
#------------------------------------------------------------------------------------------------
## Import modules
import arceditor
import os, arcpy, datetime, sys,time, os.path, glob, platform, traceback
from arcpy import env
from datetime import date, timedelta, datetime
from time import ctime
from os import listdir
from os.path import isfile, join

## Set date time variables
datetime120 = time.strftime("%Y%m%d_%H%M%S")
arcpy.env.overwriteOutput = True

## Input parameters
rasterdir = arcpy.GetParameterAsText(0)

def main():
    InFolder = rasterdir
    #Dest=arcpy.GetParameterAsText(1)
    
    arcpy.env.workspace=InFolder
    #The raster datasets in the input workspace
    in_raster_datasets = arcpy.ListRasters()
    
    # User Desktop
    userhome = os.path.expanduser('~')
    desktop = userhome+"/Desktop/"
    # Create output dir
    outdir = desktop+"RasterExtents"
    if not os.path.exists(outdir):
        os.makedirs(outdir)
    # Create output gdb, fc, and fields
    gdbname = "RasterExtents"+datetime120
    outgdb = arcpy.CreateFileGDB_management (outdir,gdbname)
    print "Output GDB location and name: "+str(outgdb)
    fcout = "ExtentsOutput_"+datetime120
    arcpy.CreateFeatureclass_management (outgdb,fcout,"POLYGON")
    gdbfc = str(outgdb)+'/'+fcout
    arcpy.AddField_management (gdbfc,"RasterName","String","","",100)
    arcpy.AddField_management (gdbfc,"RasterPath","String","","",250)
    arcpy.AddField_management (gdbfc,"CompressionType","String","","",50)
    arcpy.AddField_management (gdbfc,"CellHeight","Double","38","8","")
    arcpy.AddField_management (gdbfc,"CellWidth","Double","38","8","")
    arcpy.AddField_management (gdbfc,"Format","String","","",50)
    arcpy.AddField_management (gdbfc,"BandCount","Long","","","")
    arcpy.AddField_management (gdbfc,"NoDataValue","Double","38","8","")
    arcpy.AddField_management (gdbfc,"PixelType","String","","","50")
    arcpy.AddField_management (gdbfc,"FileSize","Double","38","8","")
    print "Creating polygons and populating attributes..."
    cursor = arcpy.InsertCursor(gdbfc)
    point = arcpy.Point()
    array = arcpy.Array()
    corners = ["lowerLeft", "lowerRight", "upperRight", "upperLeft"]
    for Ras in in_raster_datasets:
        feat = cursor.newRow()  
        r = arcpy.Raster(Ras)
        for corner in corners:    
            point.X = getattr(r.extent, "%s" % corner).X
            point.Y = getattr(r.extent, "%s" % corner).Y
            array.add(point)
        array.add(array.getObject(0))
        polygon = arcpy.Polygon(array)
        feat.shape = polygon
        feat.setValue("RasterName", Ras)
        feat.setValue("RasterPath", InFolder)
        feat.setValue("CompressionType",r.compressionType)
        feat.setValue("CellHeight",r.meanCellHeight)
        feat.setValue("CellWidth",r.meanCellWidth)
        feat.setValue("Format",r.format)
        feat.setValue("BandCount",r.bandCount)
        feat.setValue("NoDataValue",r.noDataValue)
        feat.setValue("PixelType",r.pixelType)
        feat.setValue("FileSize",r.uncompressedSize)
        cursor.insertRow(feat)
        array.removeAll()
    del feat
    del cursor    
  
if __name__ == '__main__':
    try:
        main()
    except Exception, e:
        import traceback
        map(arcpy.AddError, traceback.format_exc().split("\n"))
        arcpy.AddError(str(e))
exit()

Comments

Post a Comment

Popular posts from this blog

Updating GDB_Items XML Definitions

First, run a query to identify the feature class "PhysicalName" or UID from GDB_ITEMS. Note: The physical name appears as DBNAME.SCHEMA.FCNAME. Use one of the two for the WHERE clause in the other queries. First Query USE [DBName] SELECT ITEMS.ObjectID ,ITEMS.UUID ,ITEMTYPES.Name as ItemType ,ITEMS.[PhysicalName] ,ITEMS.[Definition] FROM [sde].[GDB_ITEMS] as ITEMS INNER JOIN [sde].[GDB_ITEMTYPES] ITEMTYPES ON [ITEMS].[ Type ] = [ITEMTYPES].[UUID] WHERE [PhysicalName] like '%FC Name%' In the results, click the XML link under the Definition field. You should see something similar to the XML below. Notice the highlighted line of code. This is the "Remarks" field. The node "IsNullable" is set as "false". I want the Remarks field to allow NULL values and would like to update it from "false" to "true". Note: The base and delta tables should be updated as well. XML Definition View <DEFeatu
Read more

Animation with ArcGIS Pro

Image
I collected data from Harris County Flood Control District and USGS for rainfall and stream heights.  I requested data between August 25th and September 5th, 2017 which was the time frame Hurricane Harvey rolled through the area. This animation shows the rainfall accumulation and stream height swell. The streams turn red once the height of the stream went above the top of bank. Here is the full video.
Read more

Add Field to a SQL Geodatabase using T-SQL

Image
Sometimes a new field will need to be added quickly. The usual process of testing changes and pushing said changes to production is not very time sensitive. Here is a method I use when a new field is needed ASAP. In this case the feature classes have been registered in the geodatabase and archiving has been enabled. This means there is an SDE base table and several SDE delta tables. The information needed to proceed is listed below. This will require DBO level access and SQL Server Management Studio or other suitable RDBMS interface like Toad. Database Name = DemoData New Field Name = NewField New Field Type (Length) = string(50) / nvarchar(50) SDE Field Type Code =   Unknown New Field Name = NewField Feature Class Name = Marker Base Table Name = Marker Base Table ID = Unknown Delta Table (Adds) = Unknown Delta Table (Archive) = Unknown First thing to do is find the unknown values. The SDE Field Type Code via Esri tells us that the SDE Field Type Code for String i
Read more