Python/pygrass

pygrass is a library that extends the GRASS capabilities to allow users to access the low-level GRASS API.

Workshops

XIV Italian meeting of GRASS and GFOSS users

Two workshops were presented during the XIV Italian meeting of GRASS at Genova. The workshops use ipython notebook to show the python and pygrass API, all the material are available on github (python, pygrass). All the execute examples are reported here:

• python
  • Introduction
  • Data Types
  • Syntax
  • Function
  • Class
  • Other

• pygrass
  • Modules
  • GIS
  • Vector
  • Raster

How to test library

pygrass has doctest inside its code. You can run doctest inside the North Carolina basic location, using the mapset user1.

To test the single module (file) you have to move to the source code of pygrass and run the following code (this is an example of functions.py):

 python -m doctest functions.py

Already tested modules

Working

The files are found in GRASS GIS 7, lib/python/pygrass/

• functions.py
• gis/region.py
• gis/__init__.py
• modules/__init.py
• vector/abstract.py
• vector/__init__.py
• vector/basic.py
• vector/table.py
• vector/geometry.py
• raster/abstract.py

Not working

• ...

References

• Zambelli, P., Gebbert, S., Ciolli, M., 2013. Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS). ISPRS International Journal of Geo-Information 2, 201–219. (DOI | PDF)

• pyGrass documentation

Some script for Landsat 7

#!/usr/bin/python
#Date: 7th February, 2013
#Public domain, GRASS GIS

#Run this script in the terminal as:
#python python-grass.py
#For debugging, run as:
#python -i python-grass.py
#PURPOSE
#This script processes LANDSAT 7 ETM+ images
#1 - unzip *.gz files
#2 - import files in GRASS GIS Location of your choice (r.in.gdal)
#3 - DN to Top of Atmosphere reflectance (i.landsat.toar)
#4 - TOA reflectance to Surface reflectance (i.atcorr)
#5 - NDVI (i.vi), Albedo (i.albedo), Emissivity (i.emissivity)

#USER HAS TO SET THOSE
#QUIET REPORTS
QIET=True
#OVERWRITE EXISTING FILES
OVR=False
#Define Landsat 7 sensor for i.landsat.toar
LSENSOR="tm7"
#Setup the path to the Landsat 7 Directories
rsdatapath="/home/yann/RSDATA/Myanmar/L7"
#Setup your GRASS GIS working directory
gisdb="/home/yann/GRASSDATA"
location="L7_Myanmar"
print location
mapset="PERMANENT"
#DEM input to atmospheric correction
inDEM=rsdatapath+"/dem.tif"
#set L7 Metadata wildcards
wldc_mtl="*_MTL.txt"
#wldc_met="*.met"
#Visibility distance [Km]
vis=18

#Set python path to enable finding of grass.script lib

#END OF USER CHANGES
###DO NOT CHANGE ANYTHING AFTER THIS LINE !
###
#Load necessary libraries
import os, glob, time, re
from grass import script as g
from grass.script import setup as gsetup
gisbase=os.environ['GISBASE']
gsetup.init(gisbase,gisdb,location,mapset)
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules.shortcuts import display as d
#Needed for floor()
from math import *
#env = os.environ.copy()
#env['GRASS_MESSAGE_FORMAT'] = 'gui'
#Function to get a list of L7 Directories in the rsdatapath
def fn(path):
	for top, dirs, files in os.walk(path):
		return [os.path.join(top, dir) for dir in dirs]

#START PROCESS
### PART 0: PRE-PROCESSING STUFF ###
#import DEM for atmospheric correction
#r.in.gdal(input=inDEM,output="dem",overwrite=OVR)
#r.mapcalc(expression="dem=25",overwrite=OVR)
#create a visibility map
r.mapcalc(expression="vis=18",overwrite=OVR)
#Find the central location of the Landsat file from metadata
metadata=[]
fileList=[]
L7Dirs=fn(rsdatapath)
for L7Dir in L7Dirs:
	#Ungzip all of your Landsat7 images in all your directories
#	print "Ungzip Landsat files in\t",L7Dir
#	p=os.system("gzip -d -q "+L7Dir+"/*.gz")
	#Using pthreads on multi-core machines
	#p=os.system("pigz -d "+L7Dir+"/*.gz")
	#Wait ten seconds for gzip to create the tif images
#	time.sleep(10)
	print "Import in GRASS GIS"
	for L7f in glob.glob(os.path.join(L7Dir,"*.[tT][iI][fF]")):
		f1=L7f.replace(L7Dir+"/","")
		f2=f1.replace(".TIF","")
		f3=f2.replace("_B10",".1")
		f4=f3.replace("_B20",".2")
		f5=f4.replace("_B30",".3")
		f6=f5.replace("_B40",".4")
		f7=f6.replace("_B50",".5")
		f8=f7.replace("_B61",".61")
		f9=f8.replace("_B62",".62")
		f10=f9.replace("_B70",".7")
		f11=f10.replace("_B80",".8")
		f12=f11.replace("L72","L71")
		L7r=f12.replace("_VCID_","")
		print "\t> ",L7r
		r.in_gdal(input=L7f,output=L7r,flags="e",overwrite=OVR)
		fileList.append(L7r)

	#reproject the DEM World map for the new PERMANENT location extents
	r.proj(input="dem",location="Myanmar",memory=10000,resolution=90.0,overwrite=OVR)
	#Get list of metadata files
	for metaf in glob.glob(os.path.join(L7Dir,wldc_mtl)):
		metadata.append(metaf)
	print "Metadata in:\n",metadata[0]
	with open(metadata[0],"r") as f:
		data=f.read()

	f.close()
	dt=data.split("\n")
	for idx in range(len(dt)):
		found=dt[idx].find("CORNER_UL_LON_PRODUCT = ")
		if found > 0:
			ulx=float(dt[idx].replace("CORNER_UL_LON_PRODUCT = ",""))
		found=dt[idx].find("CORNER_UL_LAT_PRODUCT = ")
		if found > 0:
			uly=float(dt[idx].replace("CORNER_UL_LAT_PRODUCT = ",""))
		found=dt[idx].find("CORNER_LR_LON_PRODUCT = ")
		if found > 0:
			lrx=float(dt[idx].replace("CORNER_LR_LON_PRODUCT = ",""))
		found=dt[idx].find("CORNER_LR_LAT_PRODUCT = ")
		if found > 0:
			lry=float(dt[idx].replace("CORNER_LR_LAT_PRODUCT = ",""))
		#Two cases, because a string starting by 0 is not converted well
		found=dt[idx].find("SCENE_CENTER_TIME = ")
		if found > 0:
			tim=dt[idx].replace("SCENE_CENTER_TIME = ","")
			print "timestamp=",str(tim)
			gmth=float(tim.split(":")[0])
			gmtm=float(tim.split(":")[1])
			gmtdec=float(tim.split(":")[2].replace("Z",""))
			print gmth,gmtm,gmtdec
		found=dt[idx].find("DATE_ACQUIRED = ")
		if found > 0:
			dat=dt[idx].replace("DATE_ACQUIRED = ","")
			date=dat.split("-")
		found=dt[idx].find("SUN_ELEVATION = ")
                if found > 0:
                        sunza=90-float(dt[idx].replace("SUN_ELEVATION = ",""))


	#L7 DN->Rad->Ref
	print "Convert DN to Rad to TOARef"
	f1=sorted(fileList)[0]
	f2=f1.replace(L7Dir+"/","")
	pref=f2[:-1]
	outpref=pref[:-2]+".toar."
	print "Prefix IN:\t",pref
	print "Prefix OUT:\t", outpref
	i.landsat_toar(input_prefix=pref,output_prefix=outpref,metfile=metadata[0],sensor=LSENSOR,quiet=QIET,overwrite=OVR)
	#Atmospheric Correction
	print "Atmospheric Correction"
	# Basic script for i.atcorr for L 7 ETM+
	#Geometrical conditions (L7ETM+)
	geom=7
	#Sensor height (satellite is -1000)
	sens_height=-1000
	#Here we suppose you have altitude (DEM) and Visibility (VIS) maps ready
	#---------------------------------------------
	#Visibility dummy value (overwritten by VIS raster input)
	vis=15
	#Altitude dummy value (in Km should be negative in this param file)
	#(overwritten by DEM raster input)
	alt=-1.200
	#datetime of satellite overpass (month, day, GMT decimal hour)
	mdh=str(int(date[1]))+" "+str(int(date[2]))+" "+str("%.2f" % (gmth+gmtm/60.0+gmtdec/3600.0))
	print "MM DD hh.ddd:\t",mdh
	# Central Lat/Long
	Long=ulx+(lrx-ulx)/2.0
	Lat=lry+(uly-lry)/2.0
	print "Center:\t(",Long, ",", Lat,")"
	#Atmospheric mode
	atm_mode=1 #Tropical
	#Aerosol model
	aerosol_mode=2 #Sri Lanka is a small island (maritime)
	#satellite band number (L5TM [25,26,27,28,29,30], L7ETM+ [61,62,63,64,65,66,67])
	satbandno=61 #Band 1 of L7ETM+ is first to undergo atmospheric correction
	#make time stamp for use in time series analysis
	dat1=str(date[2])+"-"+str(date[1])+"-"+str(date[0])
	dat=dat1.replace(" ","")
	timestamp=time.strftime("%d %b %Y",time.strptime(dat,"%d-%m-%Y"))
	print "Timestamp:\t",timestamp
	for idx in [1,2,3,4,5,7]:
		b=pref[:-2]+".toar."+str(idx)
		b_out=b.replace(".toar.",".surf.")
		param=[]
		param.append(str(geom)+" - geometrical conditions=Landsat 7 ETM+\n")
		param.append(str(mdh)+" "+str("%.2f" % Long)+" "+str("%.2f" % Lat)+" - month day hh.ddd longitude latitude (hh.ddd is in decimal hours GMT)\n")
		param.append(str(atm_mode)+" - atmospheric mode=tropical\n")
		param.append(str(aerosol_mode)+" - aerosols model=maritime\n")
		param.append(str(vis)+" - visibility [km] (aerosol model concentration), not used as there is raster input\n")
		param.append(str(alt)+" - mean target elevation above sea level [km] (here 600m asl), not used as there is raster input\n")
		param.append(str(sens_height)+" - sensor height (here, sensor on board a satellite)\n")
		param.append(str(satbandno)+" - i th band of ETM+ Landsat 7 (atcorr internal no)\n")
		f=open(os.path.join(L7Dir,"param_L7.txt"),"w")
		f.writelines(param)
		f.close()
		prm=os.path.join(L7Dir,"param_L7.txt")
		print "\t> ",b
		print "\t> ",b_out
		i.atcorr(input=b, elevation="dem", visibility="vis", parameters=prm, output=b_out, flags="ra", range=[0,1],quiet=QIET,overwrite=OVR)
		r.timestamp(map=b_out,date=timestamp,finish_=False)
		satbandno = satbandno + 1

	#Allocate surface reflectance names
	b1=pref[:-2]+".surf.1"
	b2=pref[:-2]+".surf.2"
	b3=pref[:-2]+".surf.3"
	b4=pref[:-2]+".surf.4"
	b5=pref[:-2]+".surf.5"
	b61=pref[:-2]+".toar.61"
	b62=pref[:-2]+".toar.62"
	b7=pref[:-2]+".surf.7"
	b8=pref[:-2]+".surf.8"

	### PART 1: BASIC STUFF ###
	b_in=pref[:-2]+".toar."
	b_clouds=pref[:-2]+".toar.acca"
	print "Clouds:\t\t>",b_clouds
	i.landsat_acca(input_prefix=b_in,output=b_clouds,overwrite=OVR)
	#png_clouds=L7Dir+"/"+pref[:-2]+".clouds.png"
	#d.mon(start='png',output=png_clouds,width=800,height=800)
	print "MASK:\t\tON"
	#Should always be rewritten!
	r.mask(raster=b_clouds,flags="i",overwrite=True)

	b_ndvi=pref[:-2]+".surf.ndvi"
	print "NDVI:\t",b_ndvi
	i.vi(red=b3,nir=b4,output=b_ndvi,viname="ndvi",quiet=QIET,overwrite=OVR,finish_=False)

	b_in=[b1,b2,b3,b4,b5,b7]
	b_albedo=pref[:-2]+".surf.albedo"
	print "Albedo:\t",b_albedo
	i.albedo(input=b_in,output=b_albedo,flags="lc",quiet=QIET,overwrite=OVR,finish_=False)
	
	b_emissivity=pref[:-2]+".surf.emissivity"
	print "Emissivity:\t",b_emissivity
	i.emissivity(input=b_ndvi, output=b_emissivity,quiet=QIET,overwrite=OVR,finish_=False)