[GRASS-user] Converting raster to polygon creates too many polygons
Hi, I've created subcatchments using r.watershed and I would like to convert these to polygons. My result at the moment contains too many polygons - 194. The original raster had only 9 areas, which should all have been contiguous seeing that they came from r.watershed. Here is the sequence of commands that I use, starting with r.report for the information about the subcatchments raster: - | Type of Map: raster Number of Categories: 18 | Data Type:CELL | Rows: 33443 | Columns: 33008 | Total Cells: 1103886544 |Projection: Transverse Mercator |N: -319.69982072S: -3155672.30728643 Res: 0.8826 |E: -33474.90929218W: -66483.27873573 Res: 1.1119 | Range of data:min = 2 max = 18 - # Set the resolution to 25m g.region rast=c83a_border res=25 # Convert the subcatchment raster to vector r.to.vect --overwrite input=c83a_dem_25m_clipped_basins_10_clip...@permanent output=c83a_subcatchme...@permanent feature=area v.report for the resulting vector gives me 194 categories. The original raster had only 9 areas, which should all have been contiguous seeing that they come from r.watershed. What am I doing wrong? I want an output vector with 9 polygons corresponding to the original 9 areas in the subcatchment raster. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Converting raster to polygon creates too many polygons
You might consider using -s flag to smooth the corners. Otherwise, I can't see why it doesn't work. R. Sab Hi, I've created subcatchments using r.watershed and I would like to convert these to polygons. My result at the moment contains too many polygons - 194. The original raster had only 9 areas, which should all have been contiguous seeing that they came from r.watershed. Here is the sequence of commands that I use, starting with r.report for the information about the subcatchments raster: - | Type of Map: raster Number of Categories: 18 | Data Type:CELL | Rows: 33443 | Columns: 33008 | Total Cells: 1103886544 |Projection: Transverse Mercator |N: -319.69982072S: -3155672.30728643 Res: 0.8826 |E: -33474.90929218W: -66483.27873573 Res: 1.1119 | Range of data:min = 2 max = 18 - # Set the resolution to 25m g.region rast=c83a_border res=25 # Convert the subcatchment raster to vector r.to.vect --overwrite input=c83a_dem_25m_clipped_basins_10_clip...@permanent output=c83a_subcatchme...@permanent feature=area v.report for the resulting vector gives me 194 categories. The original raster had only 9 areas, which should all have been contiguous seeing that they come from r.watershed. What am I doing wrong? I want an output vector with 9 polygons corresponding to the original 9 areas in the subcatchment raster. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Converting raster to polygon creates too many polygons
On 09/29/2010 02:11 PM, Hanlie Pretorius wrote: Hi, I've created subcatchments using r.watershed and I would like to convert these to polygons. My result at the moment contains too many polygons - 194. The original raster had only 9 areas, which should all have been contiguous seeing that they came from r.watershed. Here is the sequence of commands that I use, starting with r.report for the information about the subcatchments raster: - | Type of Map: raster Number of Categories: 18 | Data Type:CELL | Rows: 33443 | Columns: 33008 | Total Cells: 1103886544 1.1 billion cells. Do you really need that ? |Projection: Transverse Mercator |N: -319.69982072S: -3155672.30728643 Res: 0.8826 |E: -33474.90929218W: -66483.27873573 Res: 1.1119 | Range of data:min = 2 max = 18 It's not 9 areas but 16 categories. - # Set the resolution to 25m g.region rast=c83a_border res=25 I don't think changing resolution here will make any difference. (Unless you rerun r.watershed after) # Convert the subcatchment raster to vector r.to.vect --overwrite input=c83a_dem_25m_clipped_basins_10_clip...@permanent output=c83a_subcatchme...@permanent feature=area v.report for the resulting vector gives me 194 categories. The original raster had only 9 areas, which should all have been contiguous seeing that they come from r.watershed. You might be getting many small areas around the edges that are not connected with any of the subcatchments. Displaying the resulting vector should show you where the 194 areas are. You might consider getting rid of these edge areas with v.clean c83a_subcatchments tool=rmsa thresh= out=c83a_subcatchments_clean (choose a threshold large enough to catch all the edge areas, but smaller than any of the real subcatchments) What am I doing wrong? I want an output vector with 9 polygons corresponding to the original 9 areas in the subcatchment raster. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user This mail was received via Mail-SeCure System. -- Micha Silver Arava Development Co. +972-52-3665918 http://surfaces.co.il ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Converting raster to polygon creates too many polygons
Hanlie Pretoriuswrote: Hi, I've created subcatchments using r.watershed and I would like to convert these to polygons. My result at the moment contains too many polygons - 194. The original raster had only 9 areas, which should all have been contiguous seeing that they came from r.watershed. Here is the sequence of commands that I use, starting with r.report for the information about the subcatchments raster: - | Type of Map: raster Number of Categories: 18 | Data Type: CELL | Rows: 33443 | Columns: 33008 | Total Cells: 1103886544 | Projection: Transverse Mercator | N: -319.69982072 S: -3155672.30728643 Res: 0.8826 | E: -33474.90929218 W: -66483.27873573 Res: 1.1119 | Range of data: min = 2 max = 18 - In GRASS, you have full control of the computational region, consequently you must be able to explain (to the client/supervisor) every aspect of the computational region. Why is the horizontal resolution not only slightly off from exactly one meter, but also different for north-south and east-west? If these are the settings of the input DEM, why was the DEM created with these settings? If these resolution settings are different from the DEM's settings, there is again lots of explaining and justification to do. Life is much easier with, in this case, exactly 1m resolution which might need to be set when preprocessing and creating the high-res DEM. The resolution settings are just slightly different from 1m which can cause rounding errors later on (export, import, file formats storing resolution in single and not double precision, other GIS software using single and not double precision for resolution, etc.) # Set the resolution to 25m g.region rast=c83a_border res=25 Note that this gives you nsres = 24.99447494 and ewres = 25.00634049, see comment above. Also note that after converting to vector, the vector will not match the catchment raster because of the changed resolution. If you are interested in catchments at 25 meter resolution, resample the DEM to 25 meter, set the computational region now, don't touch it later, run r.watershed again, convert catchments to vector. # Convert the subcatchment raster to vector r.to.vect --overwrite input=c83a_dem_25m_clipped_basins_10_clip...@permanent output=c83a_subcatchme...@permanent feature=area If you use r.to.vect -v, the categories will correspond to cell values, making it easier to evaluate the vector. v.report for the resulting vector gives me 194 categories. The original raster had only 9 areas, which should all have been contiguous seeing that they come from r.watershed. The following 5x5 matrix, theoretical example for two basins 1 2 2 2 2 2 1 2 2 2 2 2 1 2 2 2 2 2 1 2 2 2 2 2 2 would be converted not to 2 but to 5 vector areas because the cells with value 1 are rendered as 4 separate square areas, if you use r.to.vect -v, with the same category i.e. belonging to the same basin. BTW, the pattern in the above sample matrix occurs more often with SFD than with MFD, particularly on such high resolution (roundabout 1m). Further on, I observed r.to.vect to create redundant areas which can be removed with v.dissolve. What am I doing wrong? Not much apart from messing with the resolution. g.region res=1 -a or g.region align=c83a_dem_25m_clipped_basins_10_clipped to cleanly align a clipped region to a raster map ;-) I want an output vector with 9 polygons corresponding to the original 9 areas in the subcatchment raster. You will likely get a vector map with more than 9 areas, but you can get a vector map with exactly 9 different categories. Markus M ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Converting raster to polygon creates too many polygons
2010/9/29, Markus Metz markus.metz.gisw...@googlemail.com: In GRASS, you have full control of the computational region, consequently you must be able to explain (to the client/supervisor) every aspect of the computational region. Why is the horizontal resolution not only slightly off from exactly one meter, but also different for north-south and east-west? If these are the settings of the input DEM, why was the DEM created with these settings? If these resolution settings are different from the DEM's settings, there is again lots of explaining and justification to do. Life is much easier with, in this case, exactly 1m resolution which might need to be set when preprocessing and creating the high-res DEM. The resolution settings are just slightly different from 1m which can cause rounding errors later on (export, import, file formats storing resolution in single and not double precision, other GIS software using single and not double precision for resolution, etc.) Thanks for the extensive comments. I don't really know why the resolution in the x and y directions are not the same. Even when I set them explicitly, I get a slight difference: - g.region -p vect=c83a_bor...@permanent nsres=25 ewres=25 projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -319.69982072 south: -3155672.30728643 west: -66483.27873573 east: -33474.90929218 nsres: 24.99447494 ewres: 25.00634049 rows: 1338 cols: 1320 cells: 1766160 - The original resolution of the DEM was at 25m resolution. I just assumend that the differences in x and y resolutions are a result of the projected surface being skew with respect to the cartesian axes. I say this because vectors that looked square in a latlong GCS display slightly skewed when projected to the Transverse Mercator PCS that I use in this location. My reason for running r.watershed at a much finer resolution was that I would get a smoother boundary when I then convert to vector. Obviously not a good idea. Markus M ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Converting raster to polygon creates too many polygons
Hanlie Pretorius wrote: I don't really know why the resolution in the x and y directions are not the same. Even when I set them explicitly, I get a slight difference: The -a flag is missing, see my previous post. - g.region -p vect=c83a_bor...@permanent nsres=25 ewres=25 projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -319.69982072 south: -3155672.30728643 west: -66483.27873573 east: -33474.90929218 nsres: 24.99447494 ewres: 25.00634049 rows: 1338 cols: 1320 cells: 1766160 - The original resolution of the DEM was at 25m resolution. I just assumend that the differences in x and y resolutions are a result of the projected surface being skew with respect to the cartesian axes. I say this because vectors that looked square in a latlong GCS display slightly skewed when projected to the Transverse Mercator PCS that I use in this location. Ah, ok, set the computational region in the target location with Transverse Mercator PCS and align with g.region -a to the desired resolution, then reproject, then run r.watershed and r.to.vect. My reason for running r.watershed at a much finer resolution was that I would get a smoother boundary when I then convert to vector. Obviously not a good idea. Someone still has to fix v.generalize for smoothing boundaries... Markus M ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] First time GRASS user -- how to import an img dataset?
Hello. I need to learn how to use GIS and I have chosen to do it with GRASS. My plan is to learn GIS/GRASS essentials during the next 12 months and start applying it for my research the 12 months after that. I have started reading the documentation available and playing with the sample data set. However, I would also like to start playing with a dataset for Puerto Rico (PR) and I am having problems to import it. How would I go about importing a new dataset? In particular, there is a landcover raster image in Erdas img format (prgap_landcover.img) that I would like to import from here: ftp://ftp.gap.uidaho.edu/products/Puerto_Rico/GISdata/Landcover/1_Land_cover_grid/ (1) What do I do with this img file? Do I first copy and paste it to my GIS Data Directory? If not, how do I go about telling GRASS where this file is (if I want to import it)? (2) I read about the commands to import raster Erdas/IMG files to GRASS (r.in.gdal -e in=image.img out=image). If I want to use this command, should my img file be in the GIS Data Directory first as I suggested in (1)? Where do I type this command? The GRASS Command Line window seems to open for an existing Location. (3) When do I create the location for my data? Do I do this before anything else? (4) When trying to create a location from the Location Wizard, I chose as the method Select coordinate system parameters from a list. I did this because there is another file in the above mentioned ftp folder with several parameters which seem to me the kind of information that the wizard asks (see an excerpt below). From there I chose lcc as the projection code. However, from here onwards I start having problems. For instance, I don't know whether to choose Datum with associated ellipsoid or Ellipsoid only. Also, the wizard asks me Central Parallel, First Standard Parallel, Second Standard Parallel, and I don't know how that relates to the two Standard_Parallel parameters given below and where do I get the other Parallel parameter from? I seem to be missing one. For the geodetic datum I think I have to choose nad83 and I think I will choose No transformation at the end of the wizard. Any suggestions will be very welcome. Thank you in advance. Saludos, Julio. Spatial_Reference_Information: Horizontal_Coordinate_System_Definition: Planar: Map_Projection: Map_Projection_Name: Lambert Conformal Conic Lambert_Conformal_Conic: Standard_Parallel: 18.03 Standard_Parallel: 18.43 Longitude_of_Central_Meridian: -66.43 Latitude_of_Projection_Origin: 17.83 False_Easting: 20.00 False_Northing: 20.00 Planar_Coordinate_Information: Planar_Coordinate_Encoding_Method: row and column Coordinate_Representation: Abscissa_Resolution: 15 Ordinate_Resolution: 15 Planar_Distance_Units: meters Geodetic_Model: Horizontal_Datum_Name: North American Datum of 1983 Ellipsoid_Name: Geodetic Reference System 80 Semi-major_Axis: 6378137.00 Denominator_of_Flattening_Ratio: 298.257222 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Open Streep Map
How to use openstreetmap data in GRASS GIS. In QGIS I have OSM plugin which allows to download, modify and upload of data. data can be editted offline and save changes back to OSM server. How to do these operation using GRASS GIS. sorry for cross posting -- Rashad ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Any GRASS users in Puerto Rico?
Hello. I need to learn how to use GIS and I have chosen to do it with GRASS. My plan is to learn GIS/GRASS essentials during the next 12 months and start applying it for my research the 12 months after that. I started reading the available documentation and also started playing with the sample data set and am now trying (unsuccessfully) to import data for PR into GRASS. I realize that I would benefit from discussing many GRASS related issues with other users and it would be great if I could find someone to discuss with these issues face to face. There is another person that will be learning GIS/GRASS with me, so we are already a big GRASS user wannabes in this campus. Anybody knows of other GRASS users in Puerto Rico we could contact? Thank you in advance. Saludos, Julio. (I apologize if you get this message twice as I first mistakenly posted it on a list intended for general GIS and OsGeo discussion). ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] First time GRASS user -- how to import an img dataset?
On 09/29/2010 06:12 PM, Mario Julio Barragan Arce wrote: Hello. I need to learn how to use GIS and I have chosen to do it with GRASS. My plan is to learn GIS/GRASS essentials during the next 12 months and start applying it for my Welcome, and good decision. research the 12 months after that. I have started reading the documentation available and playing with the sample data set. However, I would also like to start playing with a dataset for Puerto Rico (PR) and I am having problems to import it. How would I go about importing a new dataset? In particular, there is a landcover raster image in Erdas img format (prgap_landcover.img) that I would like to import from here: ftp://ftp.gap.uidaho.edu/products/Puerto_Rico/GISdata/Landcover/1_Land_cover_grid/ (1) What do I do with this img file? Do I first copy and paste it to my GIS Data Directory? If not, how do I go about telling GRASS where this file is (if I want to import it)? If the file is somewhere on your computer, then just enter the path to it with the 'input=' parameter to r.in.gdal. Something like: r.in.gdal in=/path/to/prgap_landcover.img out=prgap_landcover (BTW, I think you can actually even import into GRASS straight from the URL...) (2) I read about the commands to import raster Erdas/IMG files to GRASS (r.in.gdal -e in=image.img out=image). If I want to use this command, should my img file be in the GIS Data Directory first as I suggested in (1)? Where do I type this command? The GRASS Command Line window seems to open for an existing Location. (3) When do I create the location for my data? Do I do this before anything else? YES, before anything else. This step is crucial, and sometimes confuses new users. You *must* have a correct LOCATION defined at the start such that it matches the CRS of the original data. GRASS always insists that all maps in a LOCATION are projected in exactly the same CRS. From the metadata below it looks like this file is projected in the LCC projection matching the EPSG code of 2866. That's NAD83(HARN) Puerto Rico and Virgin Islands. You might try to first use the wizard to create a new location based on this EPSG code. That insures all the parameters get entered properly. Then, after starting GRASS in this location (also create a MAPSET - the GRASS way to keep maps filed in some directory which is logical for your work), try the r.in.gdal command. (4) When trying to create a location from the Location Wizard, I chose as the method Select coordinate system parameters from a list. I did this because there is another file in the above mentioned ftp folder with several parameters which seem to me the kind of information that the wizard asks (see an excerpt below). From there I chose lcc as the projection code. However, from here onwards I start having problems. For instance, I don't know whether to choose Datum with associated ellipsoid or Ellipsoid only. Also, the wizard asks me Central Parallel, First Standard Parallel, Second Standard Parallel, and I don't know how that relates to the two Standard_Parallel parameters given below and where do I get the other Parallel parameter from? I seem to be missing one. For the geodetic datum I think I have to choose nad83 and I think I will choose No transformation at the end of the wizard. Any suggestions will be very welcome. Thank you in advance. Saludos, Julio. Cheers. For future questions, let us know your OS and version of GRASS. And posting the output of commands like g.region -p, v.info, r.info will often help. Regards, Micha Spatial_Reference_Information: Horizontal_Coordinate_System_Definition: Planar: Map_Projection: Map_Projection_Name: Lambert Conformal Conic Lambert_Conformal_Conic: Standard_Parallel: 18.03 Standard_Parallel: 18.43 Longitude_of_Central_Meridian: -66.43 Latitude_of_Projection_Origin: 17.83 False_Easting: 20.00 False_Northing: 20.00 Planar_Coordinate_Information: Planar_Coordinate_Encoding_Method: row and column Coordinate_Representation: Abscissa_Resolution: 15 Ordinate_Resolution: 15 Planar_Distance_Units: meters Geodetic_Model: Horizontal_Datum_Name: North American Datum of 1983 Ellipsoid_Name: Geodetic Reference System 80 Semi-major_Axis: 6378137.00 Denominator_of_Flattening_Ratio: 298.257222 This mail was received via Mail-SeCure System. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user This mail was received via Mail-SeCure System. -- Micha Silver Arava Development Co. +972-52-3665918 http://surfaces.co.il ___ grass-user mailing list grass-user@lists.osgeo.org
[GRASS-user] Multiple arguments in a grass.command- Python
Greetings I have a script where I need to define a region based on 2 rasters output[0] and output[1]. I'm using thwe following expression: grass.run_command(g.region, rast = output[2] output[3], res= t_srx) But This is not correct. my question is how can I have both rasters without getting an error? (I have tried also with grass.run_command(g.region, rast = output[2],output[3], res= t_srx) and it didn't workl Thanks Luis ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Any GRASS users in Puerto Rico?
On Wed, Sep 29, 2010 at 6:20 PM, Mario Julio Barragan Arce mariojulio.barra...@upr.edu wrote: Hello. I need to learn how to use GIS and I have chosen to do it with GRASS. Welcome here! My plan is to learn GIS/GRASS essentials during the next 12 months and start applying it for my research the 12 months after that. I started reading the available documentation and also started playing with the sample data set and am now trying (unsuccessfully) to import data for PR into GRASS. Feel free to ask in this mailing list. I realize that I would benefit from discussing many GRASS related issues with other users and it would be great if I could find someone to discuss with these issues face to face. There is another person that will be learning GIS/GRASS with me, so we are already a big GRASS user wannabes in this campus. Anybody knows of other GRASS users in Puerto Rico we could contact? BTW: We have a user map here (voluntary registration): http://grass.osgeo.org/community/ with some registered users (sometimes the background is white when the JPL WMS temporarily fails) - perhaps you see some names there. Hope you get an answer here but only a tiny fraction of the GRASS users is subscribed to the mailing list. Cheers Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Multiple arguments in a grass.command- Python
Luis Lisboa wrote: I have a script where I need to define a region based on 2 rasters output[0] and output[1]. I'm using thwe following expression: grass.run_command(g.region, rast = output[2] output[3], res= t_srx) But This is not correct. my question is how can I have both rasters without getting an error? (I have tried also with grass.run_command(g.region, rast = output[2],output[3], res= t_srx) and it didn't workl You can pass lists or tuples for options which accept multiple values, e.g.: grass.run_command(g.region, rast = (output[2], output[3]), res= t_srx) or: grass.run_command(g.region, rast = [output[2], output[3]], res= t_srx) or: grass.run_command(g.region, rast = output[2:4], res= t_srx) [The last one requires that output is a list or tuple and not e.g. an array.] For a tuple, you need the parentheses to prevent the comma from being treated as an argument separator. -- Glynn Clements gl...@gclements.plus.com ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] MDOW relief shading in GRASS
Hi, I am revisiting a previous script I made in producing Multidirectional, oblique-weighted, shaded-relief. ArcGis has this script: http://blogs.esri.com/Support/blogs/mappingcenter/archive/2008/10/07/updated-hillshade-toolbox.aspx The commads in ARC 6.0.1 GRID by Rober Mark (USGS) are as follows: shade225 = hillshade(hawaii, 225, 30, shade, 5) shade270 = hillshade(hawaii, 270, 30, shade, 5) shade315 = hillshade(hawaii, 315, 30, shade, 5) shade360 = hillshade(hawaii, 360, 30, shade, 5) h00 = resample(hawaii, 1000) h01 = focalmean(h00) h02 = focalmean(h01) h03 = focalmean(h02) asp = aspect(h03) asp1 = con(isnull(asp), 293, asp) w225 = sqr(sin((asp1 - 225) div deg )) w270 = sqr(sin((asp1 - 270) div deg)) w315 = sqr(sin((asp1 - 315) div deg )) w360 = sqr(sin(asp1 div deg )) setcell minof temp = w225 * shade225 + w270 * shade270 + w315 * shade315 + w360 * shade360 shade4 = int(temp div 2) I tried porting in GRASS shell: g.region rast=$GIS_OPT_input -p # Compute hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun illumination angle r.shaded.relief map=$GIS_OPT_input shadedmap=shade225 altitude=30 azimuth=225 zmult=1 scale=1 r.shaded.relief map=$GIS_OPT_input shadedmap=shade270 altitude=30 azimuth=270 zmult=1 scale=1 r.shaded.relief map=$GIS_OPT_input shadedmap=shade315 altitude=30 azimuth=315 zmult=1 scale=1 r.shaded.relief map=$GIS_OPT_input shadedmap=shade360 altitude=30 azimuth=360 zmult=1 scale=1 # Resample DEM map r.neighbors in=$GIS_OPT_input out=res1 method=average size=$GIS_OPT_window # compute aspect map r.slope.aspect elevation=res1 format=degrees prec=float aspect=aspect zfactor=1.0 min_slp_allowed=0.0 # compute weights 225, 270, 315 and 360 r.mapcalc w225 = sin((aspect - 225)*sin(aspect - 225)) r.mapcalc w270 = sin((aspect - 270)*sin(aspect - 270)) r.mapcalc w315 = sin((aspect - 315)*sin(aspect - 315)) r.mapcalc w360 = sin((aspect)*sin(aspect)) #compute weighted r.mapcalc weightedshade = (((w225 * shade225) + (w270 * shade270) + (w315 * shade315) + (w360 * shade360))/2) r.colors map=weightedshade color=grey But I get not so good results, here's the output using nc_sp_08's elev_ned_30m (X0=MDOW, X1= default r.shaded.relief module) http://farm5.static.flickr.com/4148/5037622535_113924312b_b.jpg Anything wrong with my commands? -- cheers, maning -- Freedom is still the most radical idea of all -N.Branden wiki: http://esambale.wikispaces.com/ blog: http://epsg4253.wordpress.com/ -- ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
