[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