Explore/Manipulate Spatial Data

This has been a difficult week. As we approach the Final for GIS Programming, I am having feelings of apprehension regarding my Python coding skills. This week we worked on writing a script that would create a new geodatabase, copy data from one folder into the new geodatabase and populate a dictionary.

The reading included Chapter 6, "Exploring Spatial Data" and Chapter 7, "Manipulating Spatial Data" from our text book, by Paul A. Zandbergen.  This was all good information and appropriate sample code. Yet, I still had hurdles I could not surmount.

Creating a new geodatabase wasn't that bad:

>>>arcpy.CreateFileGDB_management("S:/GISProgramming/Module7/Results", "gcastillo.gdb") # Creates new gdb.

Copying data from one folder into the new geodatabase wasn't even that bad:

>>>arcpy.CopyFeatures_management(fc, "S:/GISProgramming/Module7/Results/gcastillo.gdb/" + fcdesc.basename)

But, I have not been able to populate the dictionary properly and then get it to print out. The "Good People" from GIS5103, Discussion Borad Module 7, even posted the following:
                           
---To populate dictionary---
county_seats.update({row[0]: row[2])
name = row[0]
population = row[2]
county_seats.update({name: population})
    print county_seats

Yet, as many times as I ran my script, this is all I could muster:

I am posting this now because I really do not like to wait til the last possible minute, and in fact, I did not procrastinate on this assignment. I just could not get the correct result. Here's to better luck next time...cheers...

Location Decisions - Where to Live

Almost there....stay on target...we are at the penultimate week for GIS4048. Well, of course, then there's the final project...but we won't discuss that right now.

For this Module we are working on how to use the Weighted Overlay Tool and conduct a thorough analysis for a couple that is moving to Alachua County.  But this is not just any couple; more on this later.

Our goal was clear; Find the perfect location for this couple's new home. So. here is where the challenge comes in. The criteria is:

The new home location must be close to North Florida Regional Medical Center (NFRMC)
The new home location must be close to UF
The new home location must be a neighborhood with a high percentage of people 40 to 49 years old and
The new home location must be a neighborhood with high home values

So, where to begin?

I put together my map beginning with a Base Map of Alachua County. After this I learned how to use the  Euclidean Distance tool and applied it to NFRMC and UF. I had to Reclassify the layers (Reclassify Tool) to make it possible to interpret and compare the results later. I also calculated the Percentage of Population Aged 40 - 49, and Homeownership, since these were also my "clients" criteria. I summarized the data and presented my map to them:

The above map could be used to discuss the four criteria with my "clients", but we are not done with this hypothetical scenario, It seems my "clients" discovered traffic in Alachua County and want to be very close to work to shorten their commute (they should see the traffic in Los Angeles- I learned to drive their). She is a doctor and he is a professor. They are smart and determined. It will be hard to please them, but I am going to try.

So, the next phase of this map creation was to add Weighted Overlays. This was really interesting to me. The Weighted Overlay Tool was very cool. This tool allows you to take layers and give more "influence" to certain data or criteria. I sat down and had a deep discussion with my "clients":

          The couple is very concerned about the commute to work. For Prof. UF, he knows he will have a set schedule that will only vary with the classes he is teaching. In fact, he expressed that when not in class, he could work from home. Mrs. Doctor, may have to go to the hospital at various times of the daydepending on medical emergencies and patient obligations. Therefore, living very close to NFRMC is the most important factor followed by homeownership. The age group was not a concern and in fact, younger or older neighbors made no difference to them by this point. Taking into consideration these factors, my weighting for my final map was:

Reclass_Own=20

reclas_Age=10

reclas_UF=10

reclas_hosp=60

I also used Model Builder to crunch the data with the tool.

In the end, I believe my "clients" will find this map very useful in determining where they will live, work and play.

Homeland Security: Protect and The Boston Marathon

We are wrapping-up our study on MEDS (Minimum Essential Datasets) and the Department of Homeland Security. This week we completed our maps using the Military Template by compiling layers of data on Critical Infrastructure around Boston for the 2012 Marathon. The overall task was to identify Critical Infrastructure, such as Hospitals, School, Dams and Airports and to analyze various Line-of-Sight locations around the finishline where security cameras could be placed. There were quite a few tasks to accomplish for this lab, such as, Analyze Data; add the MEDS data we compiled last week to our scenario map for this week; create a buffer around the event site; create a security buffer around critical infrastructure; identify and secure ingress and egress routes; generate hillshade; create surveillance points; generate viewshed; create line-of-sight graphs and a view in 3D. This last part was very cool and I learned quite a bit doing this so I wold like to list my actions and steps for this last part:

1-       To Create the line-of-sight profile graph, I had to first find where the tool was

a.       I asked myself, “Where is this?  Is this on the Draw toolbar...???” I had to Google ArcGIS Help to confirm where this “Tool” was located.

b.      I selected the Draw Arrow (next to the Drawing dropdown/pulldown window directly under the "Draw" title). I thought, “this can’t be that hard.”

c.       I thought wrong.

2-       Initially, the option was not available for me to select. That is, I could not select the “Profile Graph” on the 3D Analyst toolbar.

a.       Finally, the Create Profile Graph became available for me to select and it was not greyed out.

b.      I think I had to be in the Layout view not Data View.

3-       The Blue handles appeared when I double-clicked on the point. However, I did not need to double click. When I did this a second time for step  7 number 9. I clicked one time and the blue handles appeared and a box was placed around the line.

4-       After making my box around the surveillance point of my choice, the Graph popped up and I was able to select properties and enter a title and subtitle. I did this for several points.

5-       I exported the Graphs for later use and saved what I had before moving on to the next step.

      To create the View in 3D, I had to use ArcScene

a.       The first time I did this I did not notice that the ArcScene icon was on the 3D Analyst toolbar, so I opened ArcScene from my desktop.

7-       I added the layers as instructed and I recall thinking, I vaguely remember doing something with the Base Heights tab a long time ago...

8-      I almost missed the step, “Make sure the Factor to convert layer elevation values to scene units is set to Custom 1.0 and click OK” but I caught it just before exiting.

9-      Then it was back to ArcMap to select, copy and paste the line-of-sight from each surveillance point. This was difficult at first to get to work. However, after I did this a couple of times, all was good.

10-   This was an excellent learning point going back and forth between ArcMap and ArcScene.

a.       I saved my work and tried to export it as a 3D file (finishline_lineosight_gc.wrl). But, this did not seem to be what I wanted

b.   So, I tried again and selected export as 2D and then as a .jpg (finishline_lineosight_gc.jpg)—much better this time.

11-   After this it was time to compile my map.

This was a very involved lab that took two weeks to complete I learned a lot these past weeks and I am sure I will continue to advance in my skills and knowledge.  

Geoprocessing with Python

Are we really half-way through the semester already? Indeed time is flying. But I can definitely say we are having way too much fun. Right? This was a very cool Module. Now we are learning how to use the power of Python (not Monty-though that would be cool too) to do some ArcGIS Geoprocessing. ArcGIS and Python are closely integrated. ArcPy modules, classes, and functions give acess to all the Geoprocessing tools in ArcGIS. It was fun using some of them for this Lab.

Our task was to write a script that performs three geoprocessing functions: Add XY coordinates to the hospitals.shp layer using the AddXY tool; Create a 1000 meter buffer around the hospital features using the Buffer tool; and Dissolve the hospital buffers into a separate, single feature using the Dissolve tool. Also, after each of the three tasks above were performed, we needed to print the messages from that tool, using the GetMessages() function. Here's the run down on how I accomplished this for this week's lab:

I began by reviewing my data in ArcCatalog.

1. I copied the data from the R-drive and extracted into my Module 6 folder.
2. I opened a blank map so I could access the Arcpy window and manipulate the data I was given.
3. I set the environment (workspace) to my Module6 Data folder and prepared to begin working on my code.
4. From doing the Exercise, I surmised that to create the script I would need something like the following code:

arcpy.Buffer_analysis("hospitals.shp","S:/GISProgramming/Module6/Results/hospitals_buffer.shp","1000 METERS","","","ALL")

5. I spent time adjusting the above and researching ArcGIS Help to ensure I had this part correct.
6. I set the Local Variables for in_data & in_features and ran this part of the code.
7. Finally, I was convinved that my code would do the actions of create a buffer, and dissolve the hospital buffers into a separate, single feature.
8. Next, I needed to use the AddXY tool to insert an X and Y column at the end of my attribute table for the hospital layer. I found the format for AddXY_management.
9. Since I had set the environment to my workspace ("S:/GISProgramming/Module6/Data"), I could Execute the AddXY tool with the following code: arcpy.AddXY_management(in_features)
10. The above steps gave me a script that would perform the three geoprocessing functions: Add XY coordinates to the hospitals.shp layer; Create a 1000 meter buffer around the hospital features using the Buffer tool; and Dissolve the hospital buffers into a separate, single feature.
11. However, the Lab called for a separate line of code for the Dissolve feature, so I had to look up the Dissolve tool to use this tool separately.
12. The Dissolve tool ArcGIS Help stated that: Dissolve (Data Management) Syntax=

Dissolve_management (in_features, out_feature_class, {dissolve_field}, {statistics_fields}, {multi_part}, {unsplit_lines})

13. This was a tremendous help to me. I used the syntax from above to create my code.

The below code dissolves the hospital buffers into a separate, single feature:

fc = arcpy.Dissolve_management("S:/GISProgramming/Module6/Results/hospitals_buffer.shp","S:/GISProgramming/Module6/Results/hospitals_dissolved.shp")#, "", "","ALL")

14. By this time I had all three Geoprocessing functions working correctly, however, I did not have the code for printing the messages.
15. The last part for me: Find the code example for printing messages.
16. This felt like the longest part of this lab for me. I did my research and eventually found the example code at: http://resources.arcgis.com/en/help/main/10.2/index.html#//03q30000006s000000

GetMessage (arcpy)

import arcpy

fc = arcpy.GetParameterAsText(0)

arcpy.GetCount_management(fc)

# Print the first and last message returned by the last

#  tool executed (GetCount)

message_count = arcpy.GetMessageCount()

print(arcpy.GetMessage(0))

print(arcpy.GetMessage(message_count - 1))

17. The above was to print the first and last message of the tool. I modified the above code and after a while, I had:

# Prints Geoprocessing message

fc = arcpy.GetParameterAsText(1)

message_count = arcpy.GetMessageCount()

print(arcpy.GetMessage(0))

I added the above code after each Geoprocessing function and ran it a few times to make sure it worked properly. When I was satisfied my code was working properly, I saved my script as Mod6_gilcastillo.py and exited PythonWin. Here's my results:

Homeland Security - Prepare MEDS

This week we began our study of the Department of Homeland Security (DHS). Our task is to accomplish a Pepare MEDS exercise. MEDS stands for Minimum Essential Datasets and the idea for this was initiated by the DHS to ensure the Nation is prepared for any domestic hazard.

The National Preparedness Goal is a result of Homeland Security Presidential Directive-8 (HSPD-8) and is designed “to achieve and sustain risk-based target levels of capability to prevent, protect against, respond to, and recover from major events, and to minimize their impact on lives, property, and the economy through systematic and prioritized efforts by Federal, State, local and Tribal entities, their private and nongovernmental partners, and the general public.”  This information and more is available at:   National Preparedness Guidelines

The idea is that GIS can be very helpful in the event of a Natural or Man-made disaster and having the right types of information can save lives.  The datasets we are talking about are:

DHS Minimum Essential Data Sets (MEDS)

-Orthoimagery              -Boundaries

-Elevation                       -Structures

-Hydrography               -Land Cover

                -Transportation       -Geographic Names          

DHS is a huge Federal Agency composed of a number of  Agencies:

United States Citizenship and Immigration Services (USCIS)

United States Customs and Border Protection (CBP)

United States Coast Guard (USCG)

Federal Emergency Management Agency (FEMA)

Federal Law Enforcement TrainingCenter (FLETC)

Transportation Security Administration (TSA)

United States Immigration and Customs Enforcement (ICE)

United States Secret Service (USSS)

Directorate for Management

National Protection and Programs Directorate

Science and Technology Directorate

Domestic Nuclear Detection Office

Office of Health Affairs

Office of Intelligence and Analysis

Office of Operations Coordination

Office of Policy

An individual's safety has always depended on how well that person is prepared. It is no different when we are discussing municipalities. As we discuss GIS information, we are concerned about:  

1. The quality, accuracy and currency of all types of data about a particular location; 

2. Efficient and effective methods to optimize data sharing and interoperability between agencies and jurisdictions; 

3. Geospatial analysis to provide situational awareness at all stages of a homeland security operation.  

Therefore, DHS has promulgated guidelines that a comprehensive geospatial database be prepared, ready, available, and accessible to a community’s needs for the prevention, preparation, response, and recovery relating to any catastrophic event.  This is one of the most important components of a successful homeland security operation.

Our task is to compile a MEDS for the Boston area prior to the running of the Marathon - yes, we are traveling back in time. We are tasked to assemble the data to Prepare our map for use in a Homeland Security planning and operations analysis. I collected and re-named the layers indicated above (DHS MEDS) and ensured the projections were consistent. DHS requires that the MEDS data be in the North America Datum of 1983. 

There was a lot going on in this MEDS Prepare exercise, and while I don't want to leave anything out, I'll just hit the highlights. I started the project by preparing my environment and setting the default Geodatabase to my BostonData.gdb. This is a good way to start any project and a habit I am glad to be acquiring. After adding data to my map, I needed to Join a table to my Transportation layer. Working with Tabular data is interesting because I normally think of GIS as strictly images that comprise my map. But the tabular data can be highly useful in describing or grouping data to make the presentation more meaningful.  I took the  Census Feature Class Code (CFCC) table and Joined it to my Transportation layer to better describe the types of roads in the Boston area and then group many of those roads together to streamline my legend. Another "very cool" task here was that I used the Scale Range to hide or not display certain labels at Small Scales so that when you are "zoomed out" the map doesn't appear cluttered. Then, when you "zoom in" or select a much Larger Scale, say Street Level, the labels magically appear. I thought this was brilliant.      

Another interesting task was using the Extract by Mask Tool. This tool is found at: Spatial Analyst Tools > Extraction > Extract by Mask.  Essentially, you use one layer (extent) to extract features from another layer. I also learned that to create a color map, the input raster dataset must be a single band raster with integer values and a pixel depth of 16-bit unsigned or fewer.

One last area that was interesting to me was, converting the schema.ini file that contained all the GNIS information. The Geographic Names Information System (GNIS), was developed by the U.S. Geological Survey in corporation with the U.S. Board on Geographic Names.  GNIS data contains information about physical and cultural geographic features in the U.S. and associated areas, both current and historical (not including roads and highways). The database holds the federally recognized names of each feature and defines the location of the feature by state, county, USGS topographic map and geographic coordinates. The GNIS is the official vehicle for geographic names and it is used by the Federal Government as the source for applying geographic names to Federal maps & other printed and electronic products. The GNIS is also used to provide names data to government agencies and to the public, and provide the Geographic Names data layers to The National Map. In other words, GNIS is a very important system in the GIS world. The problem with the schema.ini file was that the information was not in rows and columns that made sense. When I opened the schema.ini file the first time I could not recognize headings, names, etc. the way it was formatted. So, I had to change the format from CSV Delimited to Delimited and add a parenthetical "pipe." This is the how I made the change in Notepad:  Format=Delimited(|).  The symbol in the parenthesis is a "pipe"; this character is found on the backslash key and to access it you use  "shift" , and then hit the backslash key.      

But that's not all folks. For some reason, after I made the change to the .ini file, my data still wasn't neatly formatted as it should have been. I tried closing and opening the file a few times and finally had to exit my remote session (disconnect from the remote server) and walk away for about 30 minutes. On a positive note, I had time to grab another cup of coffee and stretch my legs. When I returned and re-connected to the remote server, "abracadabra" the file opened in perfect form. This allowed me to use ArcCatalog, to right-click MA_FEATURES_20130404.txt and select the Create Feature Class and From XY Table to create my Geographic Names layer. 

One final note. I also took a look at the  USGS: The National Map  website to get a better idea of what is involved if I had to obtain all the MEDS information (the 7 datasets listed above) completely on my own. Wow! That would indeed be quite a time consuming task. I downloaded some of the data, just for fun and it took me quite a while to get to the data; select the appropriate data; request the data be sent to my email account; then download and extract the files on my computer. My map is not yet complete as we will continue to work on the MEDS process next week. But, here's a look at my MEDS layers as depicted in the Table of Contents (TOC). Remember, this is not yet complete, so please don't judge me.

GIS5103 Programming - Geoprocessing in ArcGIS

Here we are in week five-- time is flying fast. Learning the ins-and-outs of Python has been a challenge, but it has also been fun. This week we learned how to:


- Create a toolbox

- Create a tool using ModelBuilder

- Set Model Parameters

- Export a script from ModelBuilder

- Update a model-derived script to work outside of ArcMap

- Create a script tool

The Model, Script and Tool were all created to Clip the Soils layer to the Basin layer, then Select from the Soils layer attribute table "Not Prime Farmland." And finally, Erase the soil selection from the basin polygon. Working with ModelBuilder was quite enjoyable though exporting my Model as a Script left me with more work than I had anticipated. ModelBuilder does a great job of visually presenting you an overall picture of the task or tasks you are attempting to accomplish. 

A few notes on my process for creating my SoilErase model

1.       Before creating the model, I checked that “Overwrite the outputs of geoprocessing operations” was selected.

2.       Next, I created a toolbox

o   ArcCatalog > Toolboxes > My Toolboxes > New > Toolbox

3.       To create the model, I selected from my newly created toolbox, New> Model and named it SoilErase

o   I added the Basin and Soils shapefiles by dragging and dropping directly to the blank model.

o   I selected the Clip tool from the toolbox and dragged-dropped to the model canvas as well.

o   I used the connector tool to make the Soils shapefile the input and the Basin shapefile the clip feature

4.       Once I had the beginning of the model I ran it a few times to ensure it was going to execute.         

o   I had difficulties getting the Erase portion to work properly and I actually created three models.

o   My third model executed correctly and places new layers on my map:

§  soils_Clip.shp

§  soils_Clip_Select.shp

§  soils_Clip_Select_Erase.shp

5.  The model executed correctly as far as I could tell. Working on the python script was another matter. I had to modify the script, as outlined in the Lab instructions, each time I re-created the model. I also tried using debugging to find if I had other errors.  

So, when you export the model to Python, the code may or may not be nearly complete, so this is the real detective aspect of "Modelbuilder-to-Python code" and this is why I had to re-create my model several times. I built my Model using the Clip, Select and Erase tools, multiple times to determine where my errors were. I reset the Environment, spelled-out the file pathways, and added a line of code to defeat the "Output Feature Class already exits" error. The line of code to resolve this error is:

                                                       arcpy.env.overwriteOutput = True  

By the time I was done I could run the script in ArcGIS successfully and see the layers added to my map. When I ran the script as a tool outside of ArcGIS I was not sure of the results. That is I could not verify the new layers were added to my map even though they existed in my file folder.  Here's what my final output looked like:

Homeland Security- DC Crime Mapping

I can't believe we are on week 5! This week, we are looking at how to analyze data using the Kernel Density Tool. The subject of this week's Lab is to analyze DC Crime data by examining various categories of crime ad mapping the proximity of crimes by type and distance from various Police Stations.

I created a graph that shows all the reported crimes in the DC area:

To create this graph, my first step was to summarize the data and create a database file (dbf).
                 - Open the Attribute table, right click the column "Offense" and select Summarize
After I created the graph I experimented with different settings and selections.
                 - I decided to use the “palette” color and named the other items accordingly.
Analyzing this graph was a quick way to see how many of each type of offense was reported.
                 - I could clearly see that the top three Offenses were Theft, Theft from Auto and Burglary.

But this was only one step in the entire process to create the two maps for this week's Lab. Creating the first map, I also had the chance to use the Ring Buffer Tool. I created rings at distances of 0.5 miles, 1.0 miles and 2.0 miles. This told me how many total crimes were reported within these distances. The results were:

                                             0.5 Miles there were 668 crimes reported
                                             1.0 Miles there were 856 crimes reported
                                             2.0 Miles there were 556 crimes reported

I then used the formula:       ([Count_] / 2080) *100

To compute the percentage for each category of crime reported. 

                          At 0.5 miles the percentage of total crimes reported was: 32%

                          At 1.0 miles = 41%

                          At 2.0 miles = 27%

I looked at the data using other tools and finally created the following map:

What I learned from my analysis is that the highest crime rate is at three separate District HQs: 3D with 13% of total crime; 7D with 12% of total crime and 6D with 10%.  In fact, 6 of the top 7 crime rates are at District HQs with the only exception being the Asian Liaison unit that has 8% of the total crime. The data seems to indicate that more crime occurs in the vicinity of HQs and less crime occurs where there are Liaison units or smaller Substations.  However, more important to me is that I gained valuable experience doing this analysis using the tools I have studied in this course.

Compiling the second map was very interesting because I learned to use a knew tool: Kernel Density Tool (ArcToolbox > Spatial Analyst Tools > Density > Kernel Density). I won't try to explain the math that the tool uses; suffice it to say, this is an awesome tool and I had a blast using it. Here's the ArcGIS Help link that explains the Kernel Density Tool much better than I ever could:  

                

                                         The Awesome Kernel Density Tool

And here's my second map for this week's lab: