Wednesday, November 10, 2010

R Environment and PostgreSQL

Be The First To Comment
# Install PostgreSQL 8.4 - 32 bit and install pgJDBC Driver, OLEDB Driver from Stack Builder..........
# Install R 2.11.1 -32 bit
# Open R interface..

library(rJava)
library(DBI)
library(RJDBC)
library(RpgSQL)

.jinit(classpath=NULL,parameters=getOption("-Xmx4000m"),silent=FALSE,force.init=TRUE)
con <- dbConnect(pgSQL(), user = "username", password = "password", dbname = "db_name",host="hostname")
result<-dbGetQuery(con, 'select blah1, blah2 from "tbl_blh"')

str(result)

result

Cheers !!!

Automatically Generate KML Files from CSV file

Be The First To Comment
# Import System Modules

import sys, string, os, csv

# Enter the input and output file names here

inputFileName="PROJECT.csv"
outputFileName="PROJ.kml"

inFile=open(inputFileName)
outFile=open(outputFileName,'a')

# Enter the total number of records here
numobs = 12

header =""+"\n"+\
""+"\n"+\
""+"\n"+\
"Spatial Distribution of Students "+"\n"

# Read in the stand attributes file one line at a time

outFile.write(str(header))

for obs in range(numobs):
readline=inFile.readline()

# Select the variables
currow=readline.split(",")
student=currow[0]
year=currow[1]
rslat=currow[2]
rslong=currow[3]
advisor=currow[4]
studyarea=currow[5]

outFile.write(""+"\n"+\
""+student+""+"\n"+\
""+"Location: "+studyarea +","+"Thesis Year: "+year+","+"Thesis Advisor: "+advisor +""+"\n"+\
""+"\n"+\
""+rslong+","+rslat+"13"+""+"\n"+\
"
"+"\n"+\
"
"+"\n")
#print studentInfo
#outFile.write(str(studentInfo))

footer="
"+"\n"+\
"
"

outFile.write(str(footer))

inFile.close()
outFile.close()

Tuesday, June 1, 2010

Getting MODIS Image Automatically From FTP in Python

7 Comments
# This is a Python script that automatically downloads historical
# MODIS data from the LP DAAC FTP site
# This version should work for all of the tiled datasets
# It is currently hard-coded to downloaded specific MODIS tiles for
# the northern Great Plains & upper midwest

# Initailly historical date for Data transfer must be set on "lpdacc.txt".



import os, ftplib,sys,string


# Login information for accessing the LP DAAC FTP site
Hostname = "e4ftl01u.ecs.nasa.gov"
Username = "anonymous"
Password = "@anonymous"

# Get user inputs
# Base directory for the MODIS data

# Basedir = input("Enter the LP DAAC directory containing the dataset you want to download:")
Basedir="MOLT/MOD11A2.005"
print "The LP DAAC directory containing the dataset you want to download" +str(Basedir)


# Local directory for data storage
#Hdfdir = input("Enter the local directory where you want to store the hdf files:")
Hdfdir=r"H:\MODIS_LST_NDVI\MOD11A2\\"
print "The local directory where you want to store the hdf files" +str(Hdfdir)


# Empty lists for the collector functions
Dirlist = []
Filelist = []

mylog=open(r"H:\MODIS_LST_NDVI\MOD11A2\mylog.txt",'w',)

# Assigning Global ariables
i=0
k=0
flag=0

try:

# Define helper functions that are used to read files/subdirectories from the
# ftp site and store them as lists
def collector(line = ''):
global Dirlist
Dirlist.append(line)

def collector2(line = ''):
global Filelist
Filelist.append(line)



# Open ftp connection
ftp = ftplib.FTP(Hostname,Username,Password)

# Go to the directory containing the dataset of interest
ftp.cwd(Basedir)

# Involke the LIST ftp function, calling the collector function to store the
# results to Dirlist in list format
ftp.retrlines("LIST", collector)

# Get Directory listing only (Without including the sub directories)
mainDirlist=[]
myDirlist=[]
ftp.dir(mainDirlist.append)
myDirlist=mainDirlist[1:]
# parsing the Directory list
dirInfo=""


for mainDir in myDirlist:
#parsing the directory name only[2002.12.03]
mainDirname= mainDir[37:47]
dd=mainDirname[8:10]
mm=mainDirname[5:7]
yyyy=mainDirname[0:4]
#Extracging yyyy mm dd to compare with log file information
dirInfo=str(yyyy)+str(mm)+str(dd)
print mainDirname
print dirInfo

# Read the log file to retrive the information of latest downloaded data
logFileread=open(r"H:\MODIS_LST_NDVI\MOD11A2\lpdaac.txt",'r')
logFileread.seek(0)
logInfo=int(logFileread.read())
logFileread.close()

print"Local Drive Recent Log Dir # "+str(logInfo)
mylog.write("Local Drive Recent Log Dir # "+str(logInfo)+'\n')

# Coparing the logfile(already downloaded data) with recent datas in the ftp
if(int(logInfo) print "current path -->"+str(ftp.pwd())
if(flag==1): # if flag matches the criteria reset the counters
ftp.cwd("..")
k=0
flag=0
Filelist = []

# List all the files in the subdirectory
path=str(mainDirname)
ftp.cwd(path)
print "New path -->"+str(ftp.pwd())
# Filelist = ftp.dir()
##FTP Directory bhitra chire pni file ma chire ko chhina
ftp.retrlines("LIST", collector2)
#ftp.retrlines("LIST")
# Download data from the MODIS tiles that we are interested in
for Currow2 in Filelist:
Splitrow2 = Currow2.split()
Permissions = Splitrow2[0]

# Skip over the jpeg browse images - some of these cause problems
if Permissions[0:3] == "-rw":
Directories = Splitrow2[1]
Group = Splitrow2[2]
Size = Splitrow2[3]
Month = Splitrow2[4]
Date = Splitrow2[5]
Time = Splitrow2[6]
Filename = Splitrow2[7]
mylog.write(Filename)
LocalFile = Hdfdir + Filename
Splitfname = Filename.split(".")
# Split the header file name into its various components
Splitfname = Filename.split(".")
Mdataset = Splitfname[0]
Maqdate = Splitfname[1]
Mlocation = Splitfname[2]
Mprocdate = Splitfname[3]
Mext1 = Splitfname[4]
Mext2 = Splitfname[5]


# Pull out the horizontal and vertical tile numbers
Htile = Mlocation[1:3]
Vtile = Mlocation[4:6]


# Only retrieve data for the three tiles covering the NGP/Upper Midwest
if (((Htile == "09") & (Vtile == "04"))|((Htile == "10") & (Vtile == "04"))|((Htile == "11") & (Vtile == "04"))|((Htile == "12") & (Vtile == "04"))|((Htile == "09") & (Vtile == "05"))|((Htile == "10") & (Vtile == "05"))|((Htile == "11") & (Vtile == "05"))):
# Retrieve the hdf and xml files and place them in the local directory

ftp.retrbinary("RETR " + Filename, open(LocalFile, "wb").write)

# Write download information in the log file
logFwrite=open(r"H:\MODIS_LST_NDVI\MOD11A2\lpdaac.txt",'w')
logFwrite.seek(0)
logFwrite.write(dirInfo)
logFwrite.close()

k=k+1
if(k==14): # value of k should be no of tiles*2[for *.hdf and *.hdf.xml ]
flag=1

else:
print i
i=i+1
print "loginfor-->"+str(logInfo)
print "dirInfor-->"+str(dirInfo)
print "Dirname-->"+str(mainDirname)
print "Filename-->"+str(Filename)
mylog.write("loginfor-->"+str(logInfo)+'\n')
mylog.write("dirInfor-->"+str(dirInfo)+'\n')
mylog.write("Dirname-->"+str(mainDirname)+'\n')
mylog.write("Filename-->"+str(Filename)+'\n')



else:
print "Already downloaded in our local drive" +str(mainDirname)
mylog.write( "Already downloaded in our local drive" +str(mainDirname)+'\n')



finally:
ftp.quit()
ftp.close()
print "Closing FTP"
mylog.write("Closing FTP")
mylog.close()

Simple GUI Demo in Python

Be The First To Comment
from Tkinter import*
mainForm=Tk()

lblWelcome = Label(mainForm, text="Welcome to SSEBB Model")
lblWelcome.pack()

lblDisplay=Listbox(mainForm)

txtMyval=Entry(mainForm)
txtMyval.focus_set()
txtMyval.pack()

btnQuit=Button(mainForm,text="Close",command=quit)
btnQuit.pack(side="left",padx=10,pady=10)

def ShowVal():
lblDisplay.insert(0,txtMyval.get())


btnStart=Button(mainForm,text="Show in list",command=ShowVal)
btnStart.pack()

lblDisplay.pack(padx=20,pady=10)

mainForm.mainloop()

Friday, May 21, 2010

Parallel geoprocessing in ArcGIS Desktop.

Be The First To Comment
You will need to edit the Python script supplied to support concurrent geoprocessing of multiple processes. To make these edits you will need an understanding of how to script geoprocessing using Python.

Be aware this is just an example of an approach to parallelization of geoprocessing, there may be no benefit in the approach in your particular environment or task.

# Author: ESRI
# Date: August 2009
# Purpose: This script demonstrates a methodology for parallel geoprocessing in ArcGIS Desktop.
# The input is split into parts which are processed by separate subprocesses,
# their outputs are appended into the final result. In this specific case, the
# geoprocessing undertaken is address geocoding, but any feasible geoprocessing
# can be done by simple modification of this script.
# Look for the comment line "Begin code block for GP process of interest"

try:
import arcgisscripting, os, subprocess, sys, time, traceback
gp = arcgisscripting.create(9.3)
gp.overwriteoutput = True
startTime = time.time()
#Get the number of concurrent GP processes desired.
#It is recommended this agree with the number of your available cores.
pCount = gp.getparameter(0)

inObject = gp.getparameterastext(1)
result = gp.getcount_management(inObject)
count = int(result.getoutput(0))
gp.addmessage("Processing " + inObject + " in " + str(pCount) + " concurrent processes.")

#We will use a generic approach to splitting the input into parts, namely slicing it
#with modulo arithmetic. For example, if we want to run 4 concurrent processes
#we create 4 parts from the input using "ObjectID modulo 4 = 0 (then 1,2,3)"
#Build the splitting expressions...(Note you will need to edit this code block if using input from SDE on Oracle)
inDesc = gp.describe(inObject)
delimitedOIDName = str(gp.addfielddelimiters(inObject,str(inDesc.oidfieldname)))
wsType = str(gp.describe(str(inDesc.path)).workspacetype)
queryList = []
for q in range(pCount):
if wsType == "RemoteDatabase": #SDE
#Non-Oracle
queryList.append(delimitedOIDName + " % " + str(pCount) + " = " + str(q))
#Oracle - comment out the above line and uncomment the below line when using Oracle inputs
#queryList.append("mod(\"" + delimitedOIDName + "\"," + str(pCount) + ") = " + str(q))
elif wsType == "FileSystem": #DBase
queryList.append("mod(\"" + delimitedOIDName + "\"," + str(pCount) + ") = " + str(q))
elif wsType == "LocalDatabase" and inDesc.path.endswith(".mdb"): #Personal GDB
queryList.append(delimitedOIDName + " mod " + str(pCount) + " = " + str(q))
elif wsType == "LocalDatabase" and inDesc.path.endswith(".gdb"): #File GDB
queryList.append("mod(\"" + delimitedOIDName + "\"," + str(pCount) + ") = " + str(q))

#Create a seed name for the parts
if inDesc.name.endswith((".dbf",".shp")):
seedName = str(inDesc.name).split(".")[0]
elif inDesc.name.count(".") > 1: #Fully qualified DBMS object name
seedName = str(inDesc.name).split(".")[-1]
else:
seedName = str(inDesc.name)
#Make file GDBs at a well known path (to the system); edit this to use fast disk if you have it.
#The child scripts will write into these separate file GDBs so as to avoid lock conflicts.
#import tempfile
#wkPath = tempfile.gettempdir()
wkPath = r"C:\Temp"
gdbNameList = []
j = 0
for i in range(pCount):
gdbName = wkPath+"\parallel"+str(i)+".gdb"
while gp.exists(gdbName): #Do not clobber other parallel processes
j+=1
gdbName = wkPath+"\parallel"+str(j)+".gdb"
gdbNameList.append(gdbName)
result = gp.createfilegdb_management(wkPath,os.path.basename(gdbName))
gp.workspace = "in_memory"
#Create scripts that geoprocess the parts, write these into the file GDB directories created above
scriptPathList = []
for r in range(pCount):
gdbName = gdbNameList[r]
sName = gdbName+"/parallel.py"
scriptPathList.append(sName)
scriptLines = []
scriptLines.append("import arcgisscripting\n")
scriptLines.append("gp = arcgisscripting.create(9.3)\n")
scriptLines.append("gp.overwriteoutput = True\n")
if str(inDesc.datatype) == "Table" or str(inDesc.datatype) == "DbaseTable":
scriptLines.append("result = gp.tabletotable_conversion("+\
"\""+inObject+"\","+\
"\""+gdbName+"\","+\
"\""+seedName+str(r)+"\","+\
"\""+queryList[r]+"\")\n")
else:
scriptlines.append("result = gp.featureclasstofeatureclass_conversion("+\
"\""+inObject+"\","+\
"\""+gdbName+"\","+\
"\""+seedName+str(r)+"\","+\
"\""+queryList[r]+"\")\n")
scriptLines.append("part = result.getoutput(0)\n")
#
#Begin code block for GP process of interest:
scriptLines.append("aLocator = r\"C:\Data\Work\Product Management\Geoprocessing\TestData.gdb\NZ Locator\""+"\n")
#We could use a service. Note the \\a in \\arcgis below would be a BELL literal if not escaped as in \a (rcgis)
#scriptLines.append("aLocator = r\"GIS Servers\\arcgis on bruceh\CA_Locator.GeocodeServer\""+"\n")
scriptLines.append("fieldMap = \"Address Address;Suburb Suburb;Mailtown Mailtown;Postcode Postcode\"\n")
scriptLines.append("outFC = "+"\""+gdbName+"\\"+"Result"+str(r)+"\"\n")
scriptLines.append("result = gp.geocodeaddresses_geocoding(part,aLocator,fieldMap,outFC)\n")
scriptLines.append("result = result.getoutput(0)\n")
#End code block for GP process of interest:
#
#Make stdio aware of completion - this message will be returned to the GP stdio message pipe
scriptLines.append("print "+ "\"Finished Part" + str(r) + "\"\n")
scriptLines.append("del gp\n")
#Write script
f = open(sName,'w')
f.writelines(scriptLines)
f.flush
f.close()

#Launch the parallel processes
gp.addmessage("Launching parallel processes...\n")
processList = []
for r in range(pCount):
processList.append(subprocess.Popen([r"C:\Python25\python.exe",scriptPathList[r]],\
shell=True,\
stdout=subprocess.PIPE,\
stderr=subprocess.PIPE))
gp.addmessage("Launched process " + str(r)+ "\n")
time.sleep(2)
#
#Wait for completion
messageList = []
Done = False
while not Done:
Done = True
for p in processList:
if p.poll() is None:
Done = False
else:
messageList.append("Process " + str(processList.index(p)) + " complete at " + str(time.ctime()) + "\n")
stdOutLines = p.stdout.readlines()
for s in stdOutLines:
messageList.append(s)
stdErrorLines = p.stderr.readlines()
for e in stdErrorLines:
messageList.append(e)
if not Done: #Save 2 seconds if complete...
time.sleep(2)
for m in messageList:
gp.addmessage(m)

#Merge the results
outFC = gp.getparameterastext(2)
fm = gp.createobject("fieldmappings")
vt = gp.createobject("valuetable")
for g in gdbNameList:
gp.workspace = g
if len(gp.listfeatureclasses("Result*")) > 0:
fcDesc = gp.describe(gp.listfeatureclasses("Result*")[0])
fm.addtable(fcDesc.catalogpath)
vt.addrow(fcDesc.catalogpath)
gp.workspace = "in_memory"
result = gp.merge_management(vt,outFC,fm)
#Clean up
for g in gdbNameList:
result = gp.delete_management(g)
#Done!
endTime = time.time()
gp.addmessage("Geoprocessing duration was " + str(int(endTime - startTime)) + " seconds.")
gp.addmessage("Net performance was " + \
str(int(1.0 * count / (endTime - startTime) * 3600.0)) + " records per hour.")
del gp
except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n " + \
str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"
gp.AddError(pymsg)
msgs = "GP ERRORS:\n" + gp.GetMessages(2) + "\n"
gp.AddError(msgs)

 

© 2011 GIS and Remote Sensing Tools, Tips and more .. ToS | Privacy Policy | Sitemap

About Me