[GRASS-SVN] r53485 - in grass-addons/grass6/raster: . r.wind.sun

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Oct 19 03:25:54 PDT 2012


Author: annalisapg
Date: 2012-10-19 03:25:52 -0700 (Fri, 19 Oct 2012)
New Revision: 53485

Added:
   grass-addons/grass6/raster/r.wind.sun/
   grass-addons/grass6/raster/r.wind.sun/r.wind.sun
Log:
adding r.wind.sun - evaluates visual impact

Added: grass-addons/grass6/raster/r.wind.sun/r.wind.sun
===================================================================
--- grass-addons/grass6/raster/r.wind.sun/r.wind.sun	                        (rev 0)
+++ grass-addons/grass6/raster/r.wind.sun/r.wind.sun	2012-10-19 10:25:52 UTC (rev 53485)
@@ -0,0 +1,447 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE:	r.wind.sun
+# AUTHOR(S):	Annalisa Minelli & Ivan Marchesini
+# PURPOSE:	Calculates visual impact of aerogenerators and photovoltaic panels
+# COPYRIGHT:	(C) 2012 by the GRASS Development Team
+#
+#		This program is free software under the GNU General Public
+#		License (>=v2). Read the file COPYING that comes with GRASS
+#		for details.
+#
+#############################################################################
+#%Module
+#%  description: Calculates visual impact of aerogenerators and photovoltaic panels
+#%  keywords: visibility, photovoltaic, wind
+#%End 
+#%option
+#% key: dem
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Dem of the zone
+#% required: yes
+#%end
+#%option
+#% key: impact
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Impact output map
+#% required: yes
+#%end
+#%flag
+#%  key: w
+#%  description: Perform aerogenerators' analysis of visibility
+#%end
+#%flag
+#%  key: f
+#%  description: Perform photovoltaic analysis of visibility
+#%end
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Vector file of points where to place aerogenerators
+#% required : no
+#% guisection: Wind
+#%end
+#%option
+#% key: machine
+#% type: string
+#% gisprompt: old_file,file,input
+#% key_desc: name
+#% description: Ascii file of the aerogenerator model to simulate
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: high
+#% type: double
+#% description: Aerogenerator's height [m]
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: wind
+#% type: integer
+#% description: Wind direction in degree starting from North 
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: f
+#% type: double
+#% description: Maximum distance for computing visual impact [m]
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: windfarm2
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector map 2D of aerogenerators
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: windfarm3
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector map 3D of aerogenerators
+#% required: no
+#% guisection: Wind
+#%end
+#%option
+#% key: panels
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Vector map of panels
+#% required : no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: panels_height
+#% type: double
+#% description: Height (standard) of the panels [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: panels_width
+#% type: double
+#% description: Width (standard) of the panels [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: angle
+#% type: double
+#% description: Vertical inclination angle of the panel in degrees (above terrain surface, starting from horizontal)
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: orient
+#% type: double
+#% description: Orientation angle of the panel in degrees (starting from north, clockwise)
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: panels_center_height
+#% type: double
+#% description: Height of the panels' center on the terrain [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: resolution
+#% type: double
+#% description: Choose a resolution to work [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: min_dist_from_panel
+#% type: double
+#% description: Minimum distance for computing visual impact [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+#%option
+#% key: max_dist_from_panel
+#% type: double
+#% description: Maximum distance for computing visual impact [m]
+#% required: no
+#% guisection: Photovoltaic
+#%end
+
+import sys,os,re,numpy,math
+import grass.script as grass
+
+def main():
+	dem=options["dem"];
+	impact=options["impact"];
+	input=options["input"];
+	machine=options["machine"];
+	high=options["high"];
+	wind=options["wind"];
+	f=options["f"];
+	windfarm2=options["windfarm2"];
+	windfarm3=options["windfarm3"];
+	panels=options["panels"];
+	panels_height=options["panels_height"];
+	panels_width=options["panels_width"];
+	angle=options["angle"];
+	orient=options["orient"];
+	panels_center_height=options["panels_center_height"];
+	resolution=options["resolution"];
+	min_dist_from_panel=options["min_dist_from_panel"];
+	max_dist_from_panel=options["max_dist_from_panel"];
+	fotov=flags["f"];
+	eolic=flags["w"];
+	grass.run_command('g.gisenv',set='OVERWRITE=1');
+	if eolic:
+		res_dem=float(re.split('=',re.split('\n',grass.read_command('r.info',flags='s',map=dem))[0])[1]);
+		grass.write_command('r.mapcalc', stdin = '%s = float(%s)' % ('fdem',dem)); 
+		grass.run_command('v.in.ascii', flags='z', input=machine,output='line',format='standard',fs='|',skip='0',x='1',y='2',z='0',cat='0');
+		grass.run_command('g.copy',vect='line,face');
+		a=grass.read_command('v.info',flags='g',map='line');
+		t=float((re.split('=',re.split('\n',a)[4]))[1]);
+		b=float((re.split('=',re.split('\n',a)[5]))[1]);
+		e=float((re.split('=',re.split('\n',a)[2]))[1]);
+		w=float((re.split('=',re.split('\n',a)[3]))[1]);
+		att=(t - b ) - ( e - w )*0.577350269;
+		pala_att=(e - w)/1.5;
+		scal=float(high)/att;
+		wind=int(wind);
+		orien=180 - wind;
+		grass.run_command('g.remove',vect='line_model,face_model,pointD');
+		grass.run_command('v.transform',input='line',output='line_model',xshift='0.0',yshift='0.0',zshift='0.0',xscale=scal,yscale=scal,zscale=scal,zrot=orien,layer='1');
+		grass.run_command('v.transform',input='face',output='face_model',xshift='0.0',yshift='0.0',zshift='0.0',xscale=scal,yscale=scal,zscale=scal,zrot=orien,layer='1');
+		fdiam=( e - w ) * scal * 5;
+		grass.run_command('g.region',vect='face_model');
+		center=grass.read_command('g.region', flags='cg', vect='face_model');
+		east=float(re.split('=',(re.split('\n',center)[0]))[1]);
+		north=float(re.split('=',(re.split('\n',center)[1]))[1]);
+		try:
+			os.remove('txt');
+		except OSError:
+			print "the txt file does not exist yet but soon.. it will be created.";
+		print "to be continued..";
+		str1=[str(east-1000),' ',str(north+1000),' '];
+		str3=[str(east+1000),' ',str(north+1000),' '];
+		str5=[str(east+1000),' ',str(north-1000),' '];
+		str7=[str(east-1000),' ',str(north-1000),' '];
+		info=grass.read_command('v.info', flags='g', map='face_model');
+		add=abs(float(re.split('=',re.split('\n',info)[5])[1]));
+		grass.run_command('g.region',rast='fdem',res=res_dem);
+		grass.run_command('v.drape',input=input,type='point',rast='fdem',scale='1.0',method='nearest',output='pointD');
+		a=grass.read_command('v.to.db',flags='pc', map='pointD', type='point', layer='1', qlayer='1', option='coor', units='meters', columns='est,nord,z');
+		b=re.split('\n',a)[1:-1];	
+		n=int(re.split('\|',b[-1])[0]);
+		bibidi='';
+		bobidi='';
+		stringa_linee='';
+		stringa_facce='';
+		for i in b:
+			grass.run_command('g.region',vect=input);
+			est=float(re.split('\|',i)[1]);
+			nord=float(re.split('\|',i)[2]);
+			more=float(re.split('\|',i)[3]);
+			addmore=add+more;
+			str2=[str(est-1000),' ',str(nord+1000)];
+			str4=[str(est+1000),' ',str(nord+1000)];
+			str6=[str(est+1000),' ',str(nord-1000)];
+			str8=[str(est-1000),' ',str(nord-1000)];
+			txt=open('txt','w');
+			txt.writelines(str1+str2);
+			txt.write('\n');
+			txt.writelines(str3+str4);
+			txt.write('\n');
+			txt.writelines(str5+str6);
+			txt.write('\n');
+			txt.writelines(str7+str8);
+			k=str(re.split('\|',i)[0]);
+			palo_facce='palo_facce_'+k;
+			palo_linee='palo_linee_'+k;
+			txt='txt'
+			grass.run_command('v.transform',flags='m',input='line_model',output=palo_linee,pointsfile=txt,xshift='0.0',yshift='0.0',zshift=addmore,xscale='1.0',yscale='1.0',zscale='1.0',zrot='0.0',layer='1');			
+			grass.run_command('v.transform',flags='m',input='face_model',output=palo_facce,pointsfile=txt,xshift='0.0',yshift='0.0',zshift=addmore,xscale='1.0',yscale='1.0',zscale='1.0',zrot='0.0',layer='1');
+			if k == n :
+				stringa_linee=stringa_linee+palo_linee;
+				stringa_facce=stringa_facce+palo_facce;
+			else:
+				stringa_linee=stringa_linee+palo_linee+',';
+				stringa_facce=stringa_facce+palo_facce+',';
+			grass.run_command('v.patch',input=stringa_linee,output=windfarm2);
+			grass.run_command('v.patch',input=stringa_facce,output=windfarm3);
+		grass.run_command('g.region',vect=windfarm3,n='n'+'+'+str(f),s='s'+'-'+str(f),e='e'+'+'+str(f),w='w'+'-'+str(f));
+		tdiam=( e - w ) * scal * 3;
+		sdiam=( e - w ) * scal * 7;
+		grass.run_command('v.buffer',input=input,output='buf3',type='point,line,area',layer='1',distance=tdiam,scale='1.0',tolerance='0.01');
+		grass.run_command('v.buffer',input=input,output='buf5',type='point,line,area',layer='1',distance=fdiam,scale='1.0',tolerance='0.01');
+		grass.run_command('v.buffer',input=input,output='buf7',type='point,line,area',layer='1',distance=sdiam,scale='1.0',tolerance='0.01');
+		grass.run_command('v.patch',input='buf7,buf5,buf3',output='areas');
+		grass.run_command('nviz',elevation='fdem',vector=windfarm3);
+		pala_scal= pala_att * scal;
+		h=1;
+		for i in b:
+			grass.run_command('g.region',flags='a',vect=input,res=res_dem,n='n'+'+'+str(f),s='s'+'-'+str(f),e='e'+'+'+str(f),w='w'+'-'+str(f));
+			ca=int(re.split('\|',i)[0]);
+			xcoor=float(re.split('\|',i)[1]);
+			ycoor=float(re.split('\|',i)[2]);
+			more=float(re.split('\|',i)[3]);
+			grass.write_command('r.mapcalc', stdin = '%s = x() - %s' % ('px',xcoor));
+			grass.write_command('r.mapcalc', stdin = '%s = y() - %s' % ('py',ycoor));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((%s)^2) + ((%s)^2))' % ('dist_or','px','py'));
+			m_u= float(high) + pala_scal;
+			m_l= float(high);
+			m_h= (m_u * 0.5) + (m_l * 0.5);
+			moree= more + m_h;
+			coordinate=str(xcoor)+','+str(ycoor);
+			grass.run_command('r.viewshed',input=dem,output='step_up',coordinate=coordinate,obs_elev=m_u,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_low',coordinate=coordinate,obs_elev=m_l,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_half',coordinate=coordinate,obs_elev=m_h,max_dist=f);
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s+%s-%s))^2)+((%s)^2))' % ('dist_up',more,m_u,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s+%s-%s))^2)+((%s)^2))' % ('dist_low',more,m_l,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s-%s))^2)+((%s)^2))' % ('dist_half',moree,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_up','step_half','step_up'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_low','step_low','step_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s),%s*tan(%s))' % ('half_len',dem,moree,'dist_half','gamma_low','dist_half','gamma_up'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, abs(%s+%s-%s), abs(%s+%s-%s))' % ('disl',dem,moree,more,m_l,dem,more,m_u,dem));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((%s)^2)+((%s)^2))*(%s/%s)' % ('c_1','disl','dist_or','gamma_low','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*cos(%s), %s*cos(%s))' % ('h_1',dem,moree,'c_1','gamma_low','c_1','gamma_up'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s, %s/cos(%s))' % ('r',dem,moree,'c_1','h_1','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s, %s*%s/%s)' % ('pala',dem,moree,pala_scal,pala_scal,'r','dist_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s/%s' % ('base_1','half_len','h_1','dist_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s), %s*tan(%s))' % ('base_2',dem,moree,'h_1','gamma_up','h_1','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s+%s' % ('len','base_1','base_2'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s*%s/2' % (str(h)+'area_semi',math.pi,'len','pala'));
+			m_l= float(high) - pala_scal;
+			m_h= float(high);
+			moree= more + m_h;
+			grass.run_command('r.viewshed',input=dem,output='step_up',coordinate=coordinate,obs_elev=m_u,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_low',coordinate=coordinate,obs_elev=m_l,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_half',coordinate=coordinate,obs_elev=m_h,max_dist=f);			
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s+%s-%s))^2)+((%s)^2))' % ('dist_up',more,m_u,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s+%s-%s))^2)+((%s)^2))' % ('dist_low',more,m_l,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((abs(%s-%s))^2)+((%s)^2))' % ('dist_half',moree,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_up1','step_half','step_up'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_low','step_low','step_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s),%s*tan(%s))' % ('half_len',dem,moree,'dist_half','gamma_low','dist_half','gamma_up1'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, abs(%s+%s-%s), abs(%s+%s-%s))' % ('disl',dem,moree,more,m_l,dem,more,m_u,dem));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt(((%s)^2)+((%s)^2))*(%s/%s)' % ('c_2','disl','dist_or','gamma_low','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*cos(%s), %s*cos(%s))' % ('h_2',dem,moree,'c_2','gamma_low','c_2','gamma_up1'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s/%s' % ('pala',pala_scal,'h_2','dist_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s/%s' % ('base_1','half_len','h_2','dist_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s), %s*tan(%s))' % ('base_2',dem,moree,'h_2','gamma_up1','h_2','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s+%s' % ('len','base_1','base_2'));
+			grass.write_command('r.mapcalc', stdin = '%s = 0.5*%s*%s*%s' % (str(h)+'area_elli',math.pi,'len','pala'));
+			m_l=0.0;
+			m_h= m_u - pala_scal * 2;
+			moree= more + m_h;
+			grass.run_command('r.viewshed',input=dem,output='step_up',coordinate=coordinate,obs_elev=m_u,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_low',coordinate=coordinate,obs_elev=m_l,max_dist=f);
+			grass.run_command('r.viewshed',input=dem,output='step_half',coordinate=coordinate,obs_elev=m_h,max_dist=f);
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s+%s-%s)^2+(%s)^2)' % ('dist_up',more,m_u,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s+%s-%s)^2+(%s)^2)' % ('dist_low',more,m_l,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s-%s)^2+(%s)^2)' % ('dist_half',moree,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_up2','step_half','step_up'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s-%s' % ('gamma_low2','step_low','step_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s),%s*tan(%s))' % ('half_len',dem,moree,'dist_half','gamma_low2','dist_half','gamma_up2'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, abs(%s+%s-%s), abs(%s+%s-%s))' % ('disl',dem,moree,more,m_l,dem,more,m_u,dem));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s)^2+(%s)^2)*(%s/%s)' % ('c_3','disl','dist_or','gamma_low','gamma_low'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*cos(%s), %s*cos(%s))' % ('h',dem,moree,'c_3','gamma_low2','c_3','gamma_up2'));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s+%s-%s)^2+(%s)^2)' % ('dist_centro_ellisse',moree,pala_scal,dem,'dist_or'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s/%s' % ('base_1','half_len','h','dist_half'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*tan(%s), %s*tan(%s))' % ('base_2',dem,moree,'h','gamma_up2','h','gamma_low2'));
+			grass.write_command('r.mapcalc', stdin = '%s = (%s/%s)*if(%s <= %s, 0.5*%s, 0.5*%s)' % ('h_3','dist_centro_ellisse',pala_scal,dem,moree,'base_2','base_1'));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*%s/%s' % ('pala',pala_scal,'h_3','dist_centro_ellisse'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s <= %s, %s*0.5, %s*0.5)' % ('semiasse_ell',dem,moree,'base_2','base_1'));
+			grass.write_command('r.mapcalc', stdin = '%s = (%s*%s*%s)+((%s*(3+%s/%s))*%s*0.5)' % (str(h)+'area_composta',math.pi,'semiasse_ell','pala',scal,m_h,high,m_h));
+			patch_aree=str(h)+'area_composta'+','+str(h)+'area_elli'+','+str(h)+'area_semi';
+			grass.run_command('r.patch',flags='z',input=patch_aree,output=str(h)+'tot_area');
+			grass.run_command('r.null',map=str(h)+'tot_area',null='0');
+			grass.run_command('r.patch',input='c_3,c_2,c_1',output=str(h)+'tot_distance');
+			grass.run_command('r.null',map=str(h)+'tot_distance',null='3000000');
+			if h == n :
+				bibidi=bibidi+str(h)+'tot_distance';
+				bobidi=bobidi+str(h)+'tot_area';
+			else:
+				bibidi=bibidi+str(h)+'tot_distance'+',';
+				bobidi=bobidi+str(h)+'tot_area'+',';
+			h= h + 1;
+		grass.write_command('r.mapcalc', stdin = '%s = min(%s)' % ('dist_min',bibidi));
+		grass.run_command('r.null',map='dist_min',setnull='3000000');
+		grass.write_command('r.mapcalc', stdin = '%s = 4*%s*%s^2' % ('fov',math.pi,'dist_min'));
+		grass.write_command('r.mapcalc', stdin = '%s = 0' % ('area'));
+		grass.write_command('r.mapcalc', stdin = '%s = 0' % ('what_map'));
+		grass.write_command('r.mapcalc', stdin = '%s = 0' % ('area_others'));
+		bu=[];
+		for i in b:
+			z= 0;
+			t= 1;
+			cat=int(re.split('\|',i)[0]);
+			grass.write_command('r.mapcalc', stdin = '%s = if(min(%s) == %s, %s,%s)' % ('what_map',bibidi,str(cat)+'tot_distance',cat,'what_map'));
+			for z in re.split(',',bibidi):
+				grass.run_command('r.null',map=z,setnull='3000000');
+				grass.run_command('r.null',map=z,null='0');
+				bu.append(z);
+			z= 0;
+			c=bu[:];
+			print "c vale: ",c
+			print "cat ed n valgono rispett: ",cat,n
+			for z in c:
+				del c[cat-1];
+				del c[n-1:];
+				grass.write_command('r.mapcalc', stdin = '%s = %s+(%s*%s/%s)' % ('area_others','area_others',z,str(cat)+'tot_area',str(cat)+'tot_distance'));
+				c=bu[:];
+				print "c vale: ",c
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s = %s,%s+%s,%s+%s)' % ('area','what_map',cat,'area',str(cat)+'tot_area','area','area_others'));
+		grass.write_command('r.mapcalc', stdin = '%s = %s/%s' % ('pre','area','fov'));
+		grass.run_command('v.to.rast',input=input,output='out',use='cat',type='point',layer='1',value='1',rows='4096');
+		grass.write_command('r.mapcalc', stdin = '%s = abs(if(isnull(%s),if(%s < 1,%s*100,100),0))' % (impact,'out','pre','pre'));
+		grass.run_command('r.null',map=impact,setnull='0');
+		grass.run_command('g.remove',flags='f',rast='out,pre,area,area_others,base_1,base_2,c_1,c_2,c_3,disl,dist_centro_ellisse,dist_half,dist_low,dist_or,dist_min,dist_up,fdem,fov,gamma_low,gamma_low2,gamma_up,gamma_up1,gamma_up2,h,h_1,h_2,h_3,half_len,len,pala,px,py,r,semiasse_ell,step_half,step_low,step_up,what_map',vect='areas,buf3,buf5,buf7,face,face_model,line,line_model,pointD');
+		for i in b:
+			num=str(re.split('\|',i)[0]);
+			r2remove=num+'area_composta'+','+num+'area_elli'+','+num+'area_semi'+','+num+'tot_area'+','+num+'tot_distance';
+			v2remove='palo_linee_'+num+','+'palo_facce_'+num;
+			grass.run_command('g.remove',flags='f',rast=r2remove,vect=v2remove);
+	elif fotov:
+		grass.write_command('r.mapcalc', stdin = '%s = 0' % ('n_pann_visib'));
+		grass.write_command('r.mapcalc', stdin = '%s = 0.0' % ('dist_or'));
+		grass.write_command('r.mapcalc', stdin = '%s = 0.0' % (impact));
+		n=0;
+		grass.run_command('g.region',res=resolution);
+		grass.run_command('v.drape',input=panels,type='point,centroid,line,boundary,face,kernel',rast=dem,output='panels3D',method='nearest',scale='1.0',layer='1');
+		d=grass.read_command('v.to.db',flags='p', map='panels3D', type='point,centroid', layer='1', qlayer='1', option='coor', columns='x,y');
+		e=re.split('\n',d)[1:-1];
+		for i in e:
+			n= n + 1;
+			print n;
+			x=float(re.split('\|',i)[1]);
+			y=float(re.split('\|',i)[2]);
+			z=float(re.split('\|',i)[3]);
+			cat=int(re.split('\|',i)[0]);
+			obs_elev= 1.70 + float(panels_center_height);
+			grass.run_command('g.remove',rast='los_degree_'+str(cat));
+			coordinate=str(x)+','+str(y);
+			grass.run_command('r.viewshed',input=dem,output='los_degree_'+str(cat),coordinate=coordinate,obs_elev=obs_elev,max_dist=max_dist_from_panel);
+			grass.write_command('r.mapcalc', stdin = '%s = %s/%s' % ('los_boolean','los_degree_'+str(cat),'los_degree_'+str(cat)));
+			grass.run_command('r.null',map='los_boolean',null='0');
+			grass.write_command('r.mapcalc', stdin = '%s = %s+%s' % ('n_pann_visib','n_pann_visib','los_boolean'));
+			grass.write_command('r.mapcalc', stdin = '%s = if( %s == 0 ,null(),1)' % ('MASK','los_boolean'));
+			grass.write_command('r.mapcalc', stdin = '%s = x()-%s' % ('px',x));
+			grass.write_command('r.mapcalc', stdin = '%s = y()-%s' % ('py',y));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s)^2+(%s)^2)' % ('dist_or_'+str(cat),'px','py'));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s < %s,%s,%s)' % ('dist_or_'+str(cat),'dist_or_'+str(cat),min_dist_from_panel,min_dist_from_panel,'dist_or_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = (if(%s >= 0 & %s >= 0,90-atan(%s/%s),if(%s >= 0 & %s < 0,90-atan(%s/%s),if(%s < 0 & %s < 0,270-atan(%s/%s),if(%s < 0 & %s >= 0,270-atan(%s/%s)))))) - %s' % ('azimuth_'+str(cat),'px','py','py','px','px','py','py','px','px','py','py','px','px','py','py','px',orient));
+			grass.write_command('r.mapcalc', stdin = '%s = abs(2*(%s*0.5)*cos(%s)*(%s/(%s+(%s*0.5)*sin(%s))))' % ('apparent_width_'+str(cat),panels_width,'azimuth_'+str(cat),'dist_or_'+str(cat),'dist_or_'+str(cat),panels_width,'azimuth_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = (%s+%s)-%s' % ('dist_ver_'+str(cat),z,panels_center_height,dem));
+			grass.write_command('r.mapcalc', stdin = '%s = sqrt((%s)^2+(%s)^2)' % ('dist_'+str(cat),'dist_or_'+str(cat),'dist_ver_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = if(%s < (90+%s),%s-90+%s,if(%s = (90+%s),90,%s-90-%s))' % ('angle','los_degree_'+str(cat),angle,'los_degree_'+str(cat),angle,'los_degree_'+str(cat),angle,'los_degree_'+str(cat),angle));
+			grass.write_command('r.mapcalc', stdin = '%s = 2*((%s*0.5)*sin(abs(%s)))*(%s/(%s+((%s*0.5)*cos(%s))))' % ('apparent_height_'+str(cat),panels_height,'angle','dist_'+str(cat),'dist_'+str(cat),panels_height,'angle'));
+			grass.write_command('r.mapcalc', stdin = '%s = 2*%s*%s' % ('circle',math.pi,'dist_or_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*(tan(75))' % ('b','dist_or_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = %s*(tan(60))' % ('d','dist_or_'+str(cat)));
+			grass.write_command('r.mapcalc', stdin = '%s = (%s+%s)*%s' % ('fov_'+str(cat),'b','d','circle'));
+			grass.write_command('r.mapcalc', stdin = '%s = (%s*%s/%s)*100' % ('imp_'+str(cat),'apparent_height_'+str(cat),'apparent_width_'+str(cat),'fov_'+str(cat)));
+			grass.run_command('g.remove',rast='MASK');
+			grass.run_command('r.null',map='imp_'+str(cat),null='0');
+			grass.write_command('r.mapcalc', stdin = '%s = %s+%s' % (impact,impact,'imp_'+str(cat)));
+			grass.run_command('r.colors',flags='e',map=impact,color='rainbow');
+			grass.run_command('r.univar',map=impact);
+			r2remove='apparent_height_'+str(cat)+','+'apparent_width_'+str(cat)+','+'azimuth_'+str(cat)+','+'dist_'+str(cat)+','+'dist_or_'+str(cat)+','+'dist_ver_'+str(cat)+','+'fov_'+str(cat)+','+'imp_'+str(cat)+','+'los_degree_'+str(cat);
+			grass.run_command('g.remove',flags='f',rast=r2remove);
+		grass.run_command('r.null',map=impact,setnull='0');
+		grass.run_command('g.remove',rast='angle,b,circle,d,dist_or,los_boolean,n_pann_visib,px,py',vect='panels3D');
+	else:
+		print "You must choose to perform one of the two possible analysis, pick the flag!"	
+		
+			
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())



More information about the grass-commit mailing list