Код: Выделить всё
Extract data from .op files
# ---------------------------------------------------------------------------
# ReadFiles.py
# Description: Module used for extracting data from NCDC .op files
# Uses: Extracting lines for a given date, period or year from all files in a folder
# ---------------------------------------------------------------------------
importos, re, datetime, time, arcpy
def checkDate(date):
#date is considered correct between 1920 (first year for NCDC data is 1928)
and 2019.
#Maybe regex can be improved?
data = re.search('19[2-9]|20[0-1][0-9]([0][0-9])|([1][0-2])[0-3][0-9]', date)
ifdata == None: print 'Date format incorrect' return False
return True
def GetFileList(folder, recursive):
""" Extracts list of .op files from folder and subfolders"""
try:
theList = [] Files = os.listdir(folder) foraFile inFiles:
fullPath = folder + '\\'+ aFile #Enter all subdirectory tree recursive ifos.path.isdir(fullPath) andrecursive:
theList += GetFileList(fullPath, recursive) else: #If it's an .op file, should add to list iffullPath.endswith('.op'):
theList.append(fullPath) except: print "Error in Getting File List from the folder"
arcpy.AddError("Error in Getting File List from the folder")
arcpy.AddMessage("Error: "+ str(arcpy.GetMessages(2))) return-1
returntheList
def ExtractDate(path, date):
"""Extracts line which corresponds to given date from .op File"""
"""Returns tuple with line elements"""
#Try to convert date to datetime
try: #dateT = datetime.datetime(*time.strptime(date,"%d/%m/%Y")[0:6])
- 253 -
dateT2 = datetime.datetime.strptime(date, "%d/%m/%Y")
exceptValueError: raise "Date value not correct"
#Date is ok, open file for reading and start searching for line with date
try:
opFile = open(path, 'r')
line = opFile.readline() whileline.find(dateT2.strftime("%Y%m%d")) == -1:
line = opFile.readline() ifline == '': break else: #if loop not interrupted with break (= EOF) then line was found
opFile.close() delopFile
fields = tuple(line.split()) returnfields #if line was not found return -1 and do nothing else
opFile.close() delopFile return-1
except: raise "Error in Extract Date"
def ExtractPeriod(path, dateStart, dateEnd):
"""extracts all lines for a given period to a list"""
"""!!the period must lie within the same year because the function only
processes a opFile at a time"""
#datestart and dateend in format YYYYmmdd
#create datetime objects from both dates
try: #dates are ok, iterate dates
dateSt = datetime.datetime(*time.strptime(dateStart, "%m/%d/%Y")[0:5])
dateEn = datetime.datetime(*time.strptime(dateEnd, "%m/%d/%Y")[0:5]) #verify that dateSt is before dateEn ifdateSt > dateEn: print "Dates are in incorrect order"
arcpy.AddError("Dates are in incorrect order") return-1
except: raise "Error in Date Values for Period"
try:
d = dateSt delta = datetime.timedelta(days=1)
lineList = [] #Open opFile for reading and start searching for line with start date
opFile = open(path, 'r')
line = opFile.readline() #First search for date start whileline.find(d.strftime("%Y%m%d")) == -1 andline != "":
line = opFile.readline() ifline != '': #if loop not finished at EOF then line was found #start going through days until you get to dateend whiled <= dateEn andline != '': #If I find date, I add the tuple with data from line #and pass to the next day until I get to end date or EOF
ifline.find(d.strftime("%Y%m%d")) != -1:
lineList.append(tuple(line.split())) line = opFile.readline() d += delta else: #if EOF reached and line was not found
opFile.close() delopFile return-1 #all lines read, clean after me
opFile.close() delopFile iflineList != []: returnlineList return-1
except: raise "Error in Extract Period"
def ExtractYear(path):
try:
lineList = [] #Open opFile for reading and get all lines into memory
opFile = open(path, 'r') #first line are tabel headers, we only need the values
line = opFile.readline() line = opFile.readline() whileline != '':
lineList.append(tuple(line.split())) line = opFile.readline() #all lines read, clean after me
opFile.close() delopFile iflineList != []: returnlineList return-1
except:
raise "Error in Extract Year"
Код: Выделить всё
Scripts for shapefile manipulation
# ---------------------------------------------------------------------------
# ShapefileOps.py
# Description: Module containing methods for all the needed shapefile operations
# Uses: shapefile creation, adding points for each station, assigning attributes from
a field or all fields
# ---------------------------------------------------------------------------importarcpy, re, os
importConversions
def getStationLocation(StationList, station):
"""Get location for station with the code given"""
try: #check code, 4 to 6 numbers
data = re.search('[0-9]{3,6}', str(station)) ifdata == None: print "Station format incorrect"
arcpy.AddMessage("Station format incorrect") return 0, 0 #open station lines and search for station
lines = open(StationList, 'r')
line = lines.readline() whileline.find(str(station)) == -1:
line = lines.readline() ifline == '': break else: #if loop not interrupted with break (= EOF) then line was found
lines.close() atribute = re.split(',', line)
lat = atribute[6]
lon = atribute[7] ifisinstance(lat, str) andisinstance(lon, str):
latitude = float(lat) longitude = float(lon) else: return 0, 0 iflatitude != -99.999: returnlatitude, longitude #interrupted with break (line not found #or latitude not correct) return 0, 0
except: print "Error in getlocation, Lat/Lon" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in getlocation, Lat/Lon")
arcpy.AddMessage("Error"+ arcpy.GetMessages(2)) if notlines.closed:
lines.close() return 0, 0
def CreateShapefile_Field(lines, shapefile, Field):
"""Create a shapefile with required field"""
try: print "Creating point feature class. Field: "+ Field
arcpy.AddMessage("Creating point feature class. Field: "+ Field) #first create shapefile if list has more than 0 elements iflen(lines) == 0: return-1 else: #Check if shapefile exists and delete if needed ifos.path.exists(shapefile):
os.remove(shapefile) arcpy.AddMessage("Path already exists, overwrite"+
shapefile) print "Path already exists, overwrite"+ shapefile
folder = os.path.dirname(shapefile) filename = os.path.basename(shapefile) #Spatial reference will be WGS84.
sr = arcpy.SpatialReference() sr.factoryCode = 4326
sr.create() #Create shapefile at defined location
arcpy.CreateFeatureclass_management(folder, filename, "POINT",
'', '', '', sr)
except:
arcpy.AddMessage("Error in Create Shapefile")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
#try to add fields
try:
arcpy.AddField_management(shapefile, "Stn", "long", "9") ifField == "TEMP":
arcpy.AddField_management(shapefile, "Temp", "float") elifField == "DEWP":
arcpy.AddField_management(shapefile, "Dewp", "float") elifField == "MAXTEMP":
arcpy.AddField_management(shapefile, "maxTemp", "float") elifField == "MINTEMP":
arcpy.AddField_management(shapefile, "minTemp", "float") elifField == "PRCP":
arcpy.AddField_management(shapefile, "Prcp", "float") elifField == "SNDP":
arcpy.AddField_management(shapefile, "Sndp", "float")
except: print "Error in Add Field" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in Add Field")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
def CreateShapefile(lines, shapefile):
"""Create shapefile to store points in"""
try: #first create shapefile if list has more than 0 elements iflen(lines) == 0: return-1 else: #Only create feature class if list has items print "Creating point feature class"
arcpy.AddMessage("Creating point feature class")
arcpy.AddMessage("Number of elements: "+ str(len(lines))) print "Number of elements: "+ str(len(lines))
ifos.path.exists(shapefile):
arcpy.Delete_management(shapefile) arcpy.AddMessage("Path already exists, overwrite"+
shapefile)
folder = os.path.dirname(shapefile) filename = os.path.basename(shapefile) #Spatial reference will be WGS84. Probably will allow user to
choose spatial reference?
sr = arcpy.SpatialReference() sr.factoryCode = 4326
sr.create() #Create shapefile at defined location
arcpy.CreateFeatureclass_management(folder, filename, "POINT",
'', '', '', sr)
except:
arcpy.AddMessage("Error in Create Shapefile")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
#try to add fields
try:
arcpy.AddField_management(shapefile, "Stn", "long", "9")
arcpy.AddField_management(shapefile, "Data", "long", "9")
arcpy.AddField_management(shapefile, "Temp", "float")
arcpy.AddField_management(shapefile, "Dewp", "float")
arcpy.AddField_management(shapefile, "maxTemp", "float")
arcpy.AddField_management(shapefile, "minTemp", "float")
arcpy.AddField_management(shapefile, "Prcp", "float")
arcpy.AddField_management(shapefile, "Sndp", "float")
except: print "Field already exists or could not be created" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Field already exists or could not be created")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
def getDataField(line, field):
"""Extract all data from line and return to caller"""
try: #line is a tuple so I just iterate it.
arcpy.AddMessage(str(line)) Date = int(line[2])
Stn = int(line[0]) iffield == "TEMP": ifline[3] != '9999.9': #On days with no measurements Temp remains null.
Value = Conversions.Fahrenheit_to_Celsius(float(line[3])) else:
Value = -1.0 eliffield == "DEWP": ifline[5] != '9999.9': #On days with no measurements Dewp remains null.
Value = float(line[5]) else:
Value = -1.0 eliffield == "MAXTEMP": ifline[17] != '9999.9': #On days with no measurements maxTemp remains null. ifline[17][-1:] == '*':
Value = Conversions.Fahrenheit_to_Celsius(float(line[17][:-1])) else:
Value = Conversions.Fahrenheit_to_Celsius(float(line[17])) else:
Value = -1.0 eliffield == "MINTEMP": ifline[18] != '9999.9': #On days with no measurements minTemp remains null. ifline[18][-1:] == '*':
Value = Conversions.Fahrenheit_to_Celsius(float(line[18][:-1])) else:
Value = Conversions.Fahrenheit_to_Celsius(float(line[18])) else:
Value = -1.0 eliffield == "PRCP": ifline[19] != '99.99': #On days with no measurements precipitaton remains null. ifre.match('[A-Z]', line[19][-1:]) == None:
Value = Conversions.Inch_to_mm(float(line[19])) else:
Value = Conversions.Inch_to_mm(float(line[19][:-1])) else:
Value = -1.0 eliffield == "SNDP": iffloat(line[20]) != '999.9':
Value = Conversions.Inch_to_mm(float(line[20])) else:
Value = -1.0 else:
Value = -1.0 returnDate, Stn, Value
except: print "Error in Getting Data from line" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in Getting Data from line")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
def getData(line, Stn, Date, Temp, Dewp, maxTemp, minTemp, Prcp, Sndp):
"""Extract all data from line and return to caller"""
try:
Stn = int(line[0])
Date = int(line[2]) ifline[3] != '9999.9': #On days with no measurements Temp remains null.
Temp = Conversions.Fahrenheit_to_Celsius(float(line[3])) else:
Temp = -1.0 ifline[5] != '9999.9': #On days with no measurements Dewp remains null.
Dewp = float(line[5]) else:
Dewp = -1.0 ifline[17] != '9999.9':
#On days with no measurements maxTemp remains null. ifline[17][-1:] == '*':
maxTemp = Conversions.Fahrenheit_to_Celsius(float(line[17][:-1])) else:
maxTemp = Conversions.Fahrenheit_to_Celsius(float(line[17])) else:
maxTemp = -1.0 ifline[18] != '9999.9': #On days with no measurements minTemp remains null. ifline[18][-1:] == '*':
minTemp = Conversions.Fahrenheit_to_Celsius(float(line[18][:-1])) else:
minTemp = Conversions.Fahrenheit_to_Celsius(float(line[18])) else:
minTemp = -1.0 ifline[19] != '99.99': #On days with no measurements precipitaton remains null. Any other idea? ifre.match('[A-Z]', line[19][-1:]) == None:
Prcp = Conversions.Inch_to_mm(float(line[19])) else:
Prcp = Conversions.Inch_to_mm(float(line[19][:-1])) else:
Prcp = -1.0 iffloat(line[20]) != 999.9:
Sndp = Conversions.Inch_to_mm(float(line[20])) else:
Sndp = -1.0 #return tuple with values returnStn, Date, Temp, Dewp, maxTemp, minTemp, Prcp, Sndp
except: print "Error in Getting Data from line" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in Getting Data from line")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
def AddPointsField(lines, StationList, shapefile, Field):
"""Adds points with only one field to shapefile.
Elems contains list of lines with data to be added""" try: #Insert Cursor print "Inserting points in feature class"
arcpy.AddMessage("Inserting points in feature class")
Point_Cursor = arcpy.InsertCursor(shapefile) #Create point to add to list
arcpy.AddMessage("Point Object")
Point_Object = arcpy.CreateObject('Point')
except: print "Error in CreateShapefile, insert cursor and point object" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in CreateShapefile, insert cursor and point
object")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
try:
forlinie inlines: #A line from that date. #Must find station location
station = linie[0]
lat, lon = getStationLocation(StationList, station) #if function returned 0,0 station was not found if notlat == 0.0 and notlon == 0.0: #Create and insert point in shapefile
arcpy.AddMessage("Point: "+ str(lat) + " ; "+ str(lon))
Point_Object.X = lon Point_Object.Y = lat theValue = linie[1] #if all fields have data. iftheValue != -1.0: #Create New point to add to shapefile
NewPoint = Point_Cursor.newRow() NewPoint.Shape = Point_Object #populate fields
arcpy.AddMessage("Value: "+ str(theValue))
NewPoint.Stn = long(station) ifField == "TEMP":
NewPoint.Temp = float(theValue) elifField == "DEWP":
NewPoint.Dewp = float(theValue) elifField == "MAXTEMP":
NewPoint.maxTemp = float(theValue) elifField == "MINTEMP":
NewPoint.minTemp = float(theValue) elifField == "PRCP":
NewPoint.Prcp = float(theValue) elifField == "SNDP":
NewPoint.Sndp = float(theValue) #Add point to shapefile
Point_Cursor.insertRow(NewPoint) #clean up delPoint_Cursor delPoint_Object
except: print "Error in CreateShapefile, add points" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in CreateShapefile, add points")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
def AddPoints(lines, StationList, shapefile):
"""Adds all points in lines set to shapefile""" print "Inserting points in feature class"
arcpy.AddMessage("Inserting points in feature class")
Stn, Date = 0, 0
Temp, Dewp, maxTemp, minTemp, Prcp, Sndp = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
try: #Insert Cursor
arcpy.AddMessage("Insert Cursor")
Point_Cursor = arcpy.InsertCursor(shapefile) #Create point to add to list
arcpy.AddMessage("Point Object")
Point_Object = arcpy.CreateObject('Point')
except: print "Error in CreateShapefile, insert cursor and point object" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in CreateShapefile insert cursor and point
object")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
try: forlinie inlines: #A line from that date. #Must find station location
station = linie[0]
arcpy.AddMessage("Find lat lon for "+ str(station))
lat, lon = getStationLocation(StationList, station) #if function returned 0,0 station was not found if notlat == 0.0 and notlon == 0.0: #Create and insert point in shapefile
arcpy.AddMessage("Point: "+ str(lat) + str(lon))
Point_Object.X = lon Point_Object.Y = lat Stn, Date, Temp, Dewp, maxTemp, minTemp, Prcp, Sndp = getData(linie, Stn, Date, Temp, Dewp, maxTemp, minTemp, Prcp, Sndp) #if all fields have data. ifTemp != -1.0 orDewp != -1.0 ormaxTemp != -1.0 orPrcp
!= -1.0 orminTemp != -1.0 orSndp != -1.0: #Create New point to add to shapefile
NewPoint = Point_Cursor.newRow() NewPoint.Shape = Point_Object #populate fields
arcpy.AddMessage(str(Stn) + " "+ str(Date) + " "+
str(Temp) + " "+ str(Dewp) + " "+ str(maxTemp) + " "+ str(minTemp) + " "+
str(Prcp) + " "+ str(Sndp))
NewPoint.Data = long(Date) NewPoint.Stn = long(Stn) NewPoint.Temp = float(Temp) NewPoint.Dewp = float(Dewp) NewPoint.maxTemp = float(maxTemp) NewPoint.minTemp = float(minTemp) NewPoint.Prcp = float(Prcp) NewPoint.Sndp = float(Sndp) #Add point to shapefile
Point_Cursor.insertRow(NewPoint) #clean up delPoint_Cursor delPoint_Object
except: print "Error in CreateShapefile, add points" print "Error"+ str(arcpy.GetMessages(2))
arcpy.AddMessage("Error in CreateShapefile, add points")
arcpy.AddMessage("Error"+ str(arcpy.GetMessages(2)))
