Plot Precip map

[rexmorgan]

This is working code. You will of course have to replace the directories used in this code with your own that will be unique to your particular situation. You also must have ArcGIS, and the newest version is 10 which will only work with python 2.6. There is plenty of comments in the code to try and help people understand what we area doing. I also much acknowledge the help of pyTony with his help in getting the csv file to work properly when spaces are involved and for Ene Uran's help for the initial bit of code that got going down the correct path. Thanks again to both of you. Link to related post Having said all of that here is the code to play with.

# --------------------------------------------------------------------------------------------------------
# GISClimateMapping.py
# Modified and Updated by Rex Morgan and Josh Barnwell; Original Max Min Precip Snow .py by Josh Barnwell
# Description: Creates GIS maps for Temperatures, and Rain and Snow Amounts.
# Much help was provided through daniweb.com at the following thread
# http://www.daniweb.com/software-development/python/threads/432830/importing-csv-into-array#post1857807
# --------------------------------------------------------------------------------------------------------
# User Input of Month, Day, and what Variable to be plotted

from time import *
import csv
from datetime import *
from numpy import *
now = date.today().isoformat()
print strftime("%Y-%m-%d %H:%M:%S", gmtime()) #used to evaluate the amount of time the program takes to run. 


print ' '
print '                     ARCGIS Climate Data Mapping Program'
print ' '

months = raw_input('Enter the 3 letter abbreviation of the month you want to plot: ')
names = 'Z:\\Desktop\\Calendar_Year_Rainfall_all.csv'
Rex = 'Z:\\Desktop\\COOP_rainfall.csv'
usermonth = "[" + months + "]"

#This is still yet to be enacted. 
#month1 = raw_imput('Enter the 3 letter abbreviation of the beginning month:') 
#month2 = raw_input('Enter the 3 letter abbreviation of the ending month:')

print ' '
print 'This will take about 2 to 5 minutes...'

#This section is not doing anything at the moment.  
if 'JAN' in months:
    month = 'January'
if 'FEB' in months:
    month = 'February'
if 'MAR' in months:
    month = 'March'
if 'APR' in months:
    month = 'April'
if 'MAY' in months:
    month = 'May'
if 'JUN' in months:
    month = 'June'
if 'JUL' in months:
    month = 'July'
if 'AUG' in months:
    month = 'August'
if 'SEP' in months:
    month = 'September' 
if 'OCT' in months:
    month = 'October'
if 'NOV' in months:
    month = 'November' 
if 'DEC' in months:
    month = 'December' 

# Import arcpy module
import arcpy
import arcpy.mapping

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

import sys, string, os

# Directories
tempdir = 'CSV\\'
directory = 'N:\\climate\\ClimateSheets\\' + str(tempdir)
TopoToR = 'R:\\ArcGIS\\default.gdb\\TopoToR'

# Local variables:
rr3precip_csv = 'Z:\\Desktop\\COOP.csv'
GGW_Counties = "R:\\fxc\\Western Region FXC Packet\\GGW.shp"
Var_layer = months
Var_shp = str(directory) + 'Var.shp'
Var_new = str(directory) + 'Varnew.shp'
Var_img = str(directory) + 'Var.img'
Var100_img = str(directory) + 'Var100.img'
Vari100 = str(directory) + 'Vari100'
Varcontour_shp = str(directory) + 'Varcontour.shp'
Output_stream_polyline_features = ""
Output_remaining_sink_point_features = ""
Output_diagnostic_file = ""
Output_parameter_file = ""
TopoToRasterExtent= "-108.891024 46.540403 -104.041596 49.000027"
spRef = r"Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
print 'Made it through Local variables'


import arcview
import sys, string, os, arcgisscripting
gp = arcgisscripting.create(10.0)
gp.SetProduct("ArcView") # or ArcEditor or ArcInfo
gp.overwriteoutput = 1
#arcpy.overwriteOutput = True
print 'Made it through import arcview'


# Process: Make XY Event Layer
# In the 'Make XY Evennt Layer' Rex is the file name, "Long" is the header in the file Rex, 
#   "Lat" is the header for latitude in the file Rex
#   Var_layer  is equal to months, which is equal to user input in the form of a 3 letter identifier 
arcpy.MakeXYEventLayer_management(Rex, "Long", "Lat", Var_layer, spRef, "")
print 'Make XY Event Layer Done'


# Process: Copy Features
arcpy.CopyFeatures_management(Var_layer, Var_shp, "", "0", "0", "0")
print 'Copy Features Done'


# Process: Topo to Raster
from arcpy import env
env.extent = "-108.891024 46.540403 -104.041596 49.000027" 
from arcpy.sa import *
outTTR = TopoToRaster([TopoPointElevation([[Var_shp, Var_layer]])], "0.003", "", "20", "", "", "ENFORCE", "CONTOUR", "", "", "", "", "", "")
outTTR.save(TopoToR)
print 'Topo To Raster Done'


# Process: Clip
arcpy.Clip_management(TopoToR, "-108.891024 46.540403 -104.041596 49.000027", Var_img, GGW_Counties, "", "ClippingGeometry")
print 'Clip Done'


# Process: Define Projection for Shapefile in order to plot record values correctly
arcpy.DefineProjection_management(Var_shp, spRef)
print 'Define Projection Done' # used during code testing can be eliminated in time
print strftime("%Y-%m-%d %H:%M:%S", gmtime()) # used during code testing can be eliminated in time

print 'Exporting Map' #testing code can be eliminated in time
import time
from datetime import date
now = date.today().isoformat()
today = date.today()

mxd1 = arcpy.mapping.MapDocument(r"R:\EventData\CoopPrecip1024.mxd")
mxd3 = arcpy.mapping.MapDocument(r"R:\EventData\CoopPrecip13.mxd")


if 'YTD' in months:
  print 'if statement'
  newtitle = arcpy.mapping.ListLayoutElements(mxd1, "TEXT_ELEMENT", "title")[0]
  for lyr in arcpy.mapping.ListLayers(mxd1):
   if lyr.name == "Var":
     if lyr.supports("LABELCLASSES"):
       lblClass = lyr.labelClasses[0] #this is the Default label class.
       lblClass.expression = usermonth
  newtitle.text = 'Annual COOP Precipitation Year to Date'
  arcpy.mapping.ExportToPNG(mxd1, r'R:\\Maps\\' + months + str(now) + '.png', resolution=100)
else:
  newtitle = arcpy.mapping.ListLayoutElements(mxd3, "TEXT_ELEMENT", "title")[0]
  for lyr in arcpy.mapping.ListLayers(mxd3):
   if lyr.name == "Var":
     if lyr.supports("LABELCLASSES"):
       lblClass = lyr.labelClasses[0] #this is the Default label class.
       lblClass.expression = usermonth
  newtitle.text = 'Monthly COOP Precipitation for ' + month
  arcpy.mapping.ExportToPNG(mxd3, r'R:\\Maps\\' + months + str(today) + '.png', resolution=100)


print 'Map is Complete and located in R:\Maps'

print strftime("%Y-%m-%d %H:%M:%S", gmtime()) # used to evaluate the amount of time the code took to run

Thanks again to all those that have helped. Once I have more to offer I will update this post. If you have any suggestions on incorporating the code from the linked post above don't hesitate to post your suggestions.