Re: [GRASS-user] 'Growing' a raster on the edges only (Paulo van Breugel)
Thanks Paulo, your suggestion of combining r.neigbors and r.patch worked. From: Paulo van Breugel p.vanbreu...@gmail.com Subject: Re: [GRASS-user] 'Growing' a raster on the edges only Message-ID: 516ba41a.4060...@gmail.com Content-Type: text/plain; charset=UTF-8; format=flowed Try r.neighbors function with the mean or medium neighborhood operation, or, if you are working with categorical values, use the 'mode' neighborhood operation. You can combine the result with the original layer using r.patch so only the no-data cells at the borders are filled with the newly created values You may also have a look at r.fillnulls, which fills in no-data cells using spline interpolation. Paulo ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Rainfall interpolation for a time series
Hi, I have rainfall measurements from 6 gauges that I want to interpolate to an areal value (a 'surface'), so that I can compare the the interpolated gauge values to a satellite rainfall estimate that covers a grid cell of 28kmx28km. Two of the gauges are outside, but close to the border of the grid cell. Therefore, I also need to clip the interpolated surface to the grid cell and to get the average of the surface value in this clipped surface. A search on Google revealed that v.vol.rst is probably a good choice for interpolating the rainfall because it takes elevation into account. However, for each rain gauge I have 730 values representing a daily measurement over two years. As output from GRASS, I need a text file with the interpolated rainfall values for each day in my time series. So, I was wondering if there is an 'easy' way to get my output without creating 730 vectors and rasters? Regards Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Rainfall interpolation for a time series
Hi Nick, Thanks for the reply. Did you do your comparison with a time series as well? If so, do you perhaps have the start of a script that you would be willing to share? Thanks Hanlie 2011/8/15, Nick Jachowski njachow...@gmail.com: Hi Hanlie, I've this before with 11 gauge stations, also using TRMM data. You'll probably want to write a script and use r.mask to mask the interpolation map to your grid cell, then do r.sum on your masked area. After that just divide by the area to get the average rainfall. Given the small area and limited number of gauges this shouldn't take very long to compute for 2 years of data. You should be able to fairly easily output your desired time series to a text file. Best, Nick On Mon, Aug 15, 2011 at 3:00 PM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I have rainfall measurements from 6 gauges that I want to interpolate to an areal value (a 'surface'), so that I can compare the the interpolated gauge values to a satellite rainfall estimate that covers a grid cell of 28kmx28km. Two of the gauges are outside, but close to the border of the grid cell. Therefore, I also need to clip the interpolated surface to the grid cell and to get the average of the surface value in this clipped surface. A search on Google revealed that v.vol.rst is probably a good choice for interpolating the rainfall because it takes elevation into account. However, for each rain gauge I have 730 values representing a daily measurement over two years. As output from GRASS, I need a text file with the interpolated rainfall values for each day in my time series. So, I was wondering if there is an 'easy' way to get my output without creating 730 vectors and rasters? Regards 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] DEM resolution from r.surf.contour
GRASS 6.4.1 on Windows XP. g.region output for 5m resolution: g.region -p vect=contours_20m@PERMANENT nsres=5 ewres=5 projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -3121229.69982072 south: -3156672.30728643 west: -67483.27873573 east: -32474.90929218 nsres: 4.9996625 ewres: 4.99976713 rows: 7089 cols: 7002 cells: 49637178 g.region output for 20m resolution: g.region -p vect=contours_20m@PERMANENT nsres=20 ewres=20 projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -3121229.69982072 south: -3156672.30728643 west: -67483.27873573 east: -32474.90929218 nsres: 20.00147148 ewres: 20.00478254 rows: 1772 cols: 1750 cells: 3101000 Thanks Hanlie 2011/5/6, Markus Neteler nete...@osgeo.org: On Wed, May 4, 2011 at 10:48 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I'd like to ask a question related to another question I was asking recently about converting contour vectors to rasters and creating a DEM from the raster contours using r.surf.contour. My supplied vector contour layer has 20m intervals between contour lines. When I convert this to raster, the smallest region resolution I can use is 5m, otherwise v.to.rast runs out of memory. Please post - GRASS version/OS - g.region -p output to better understand the problem. Markus When I use r.surf.contour to create the DEM from the raster, which region resolution makes sense then? I set it to 20m because of the contour interval, but I'm starting to think that's not the best option because the interpolation between two contour lines should be as fine as possible. Any thoughts? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.to.rast gives empty raster layer
Hi Moritz, I'm not sure what you mean by 'Maris' first part of the solution'. I saw only one part to Maris's suggestion (refining the region's resolution) because I said in my initial mail that I had set extents of the region to the contours_20m vector layer. Initially I did have a resolution of 20m and got an empty raster layer as a result. I share Micha's confusion about this. Although it doesn't make much sense to have a 20m cell size (my bad) I would have thought that v.to.rast would still have given me an output. Regards Hanlie 2011/5/3, Moritz Lennert mlenn...@club.worldonline.be: On 03/05/11 15:47, Micha Silver wrote: On 05/03/2011 04:34 PM, Moritz Lennert wrote: On 03/05/11 10:20, Micha Silver wrote: On 05/03/2011 09:11 AM, Hanlie Pretorius wrote: Hi Maris, Thanks, your suggestion worked. With the region's resolution set to 5m, I got some output. Hi Maris: Sorry for butting in on this, but I don't understand why the region resolution should make any difference to v.to.rast?? Would you mind to explain? From the v.to.rast manual: v.to.rast will only affect data in areas lying inside the boundaries of the current geographic region. Before running v.to.rast, the user should therefore ensure that the current geographic region is correctly set and that the region resolution is at the desired level. Thanks, Moritz What I'm trying to understand is the question about resolution. Would you expect, as Hanlie indicted, to get *no output* in the raster because of too course resolution? He claims that at 20m x 20m v.to.rast of a line feature created all NULL cells, but at 5m x 5m the lines did appear in the raster. ?? Hanlie didn't tell us whether he applied Maris' first part of the solution as well (but I guess he did): g.region vect=contours_20m I agree with you that resolution does not make the difference. However, I guess that Maris meant that having a 20m resolution seems too coarse for 20m contour lines, but even that is not clear as all depends on the topography of the region in question. If you have 100m altitude change for 50m horizontal change, then 20m is a bad choice of resolution, but if you have 100m altitude change for 1km horizontal change, it should be more than enough. So, in conclusion, the region's resolution setting can make a huge difference in v.to.rast, but in this particular case does not explain why there was no output. This is certainly linked to the region's extension. Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.to.rast gives empty raster layer
Oh, I see, thanks. r.stats input=contours_20m@CDNGI_DEM 1620-1623.294118 1639.764706-1643.058824 1659.529412-1662.823529 1679.294118-1682.588235 1699.058824-1702.352941 1718.823529-1722.117647 1738.588235-1741.882353 1758.352941-1761.647059 1778.117647-1781.411765 1797.882353-1801.176471 1817.647059-1820.941176 1837.411765-1840.705882 1857.176471-1860.470588 1876.941176-1880.235294 1900-1903.294118 1919.764706-1923.058824 1939.529412-1942.823529 1959.294118-1962.588235 1979.058824-1982.352941 1998.823529-2002.117647 2018.588235-2021.882353 2038.352941-2041.647059 2058.117647-2061.411765 2077.882353-2081.176471 2097.647059-2100.941176 2117.411765-2120.705882 2137.176471-2140.470588 2156.941176-2160.235294 2180-2183.294118 2199.764706-2203.058824 2219.529412-.823529 2239.294118-2242.588235 2259.058824-2262.352941 2278.823529-2282.117647 2298.588235-2301.882353 2318.352941-2321.647059 2338.117647-2341.411765 2357.882353-2361.176471 2377.647059-2380.941176 2397.411765-2400.705882 2417.176471-2420.470588 2436.941176-2440.235294 2456.705882-2460 * (Wed May 04 10:35:47 2011) Command finished (0 sec) For some reason I couldn't get the raster contours to display yesterday - another reason why I though the layer was empty. I adjusted the colour table an re-rendered the map with no effect. Today, however, it works and I can see it. Sorry about the confusion. 2011/5/4, Moritz Lennert mlenn...@club.worldonline.be: On 04/05/11 08:42, Hanlie Pretorius wrote: Hi Moritz, I'm not sure what you mean by 'Maris' first part of the solution'. I saw only one part to Maris's suggestion (refining the region's resolution) because I said in my initial mail that I had set extents of the region to the contours_20m vector layer. Initially I did have a resolution of 20m and got an empty raster layer as a result. I share Micha's confusion about this. Although it doesn't make much sense to have a 20m cell size (my bad) I would have thought that v.to.rast would still have given me an output. Rereading the original mail I see this: 2011/5/2, Hanlie Pretoriushanlie.pretor...@gmail.com: r.info for the result is: | Layer: contours_20m@CDNGI_DEM Date: Mon May 02 14:18:00 2011 | Mapset: CDNGI_DEM Login of Creator: hanlie | Location: SA_Lo_29E | DataBase: F:\grassdata | Title: Labels ( contours_20m@CDNGI_DEM ) | Timestamp: none | | | Type of Map: raster Number of Categories: 0 [snip] | Rows: 1772 | Columns: 1750 | Total Cells: 3101000 | Projection: Transverse Mercator | N: -3121229.69982072 S: -3156672.30728643 Res: 20.00147148 | E: -32474.90929218 W: -67483.27873573 Res: 20.00478254 | Range of data: min = 1620 max = 2460 [snip] As you can see from the number of categories (0), there is no information in this raster file. There are no category values, but there are data values, so your map actually is not empty. What does r.stats contours_20m@CDNGI_DEM give you ? 1) v.to.rast does create category labels by default. If you use the 'labelcol' parameter, then you get category values. Example in the NC demo dataset: v.to.rast elev_ned10m_cont10m out=contours col=level r.info contours: Number of Categories: 0 Range of data:min = 60.00 max = 150.00 2) v.to.rast elev_ned10m_cont10m out=contours col=level labelcol=cat r.info contours: Number of Categories: 150 Range of data:min = 60.00 max = 150.00 Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] DEM resolution from r.surf.contour
Hi, I'd like to ask a question related to another question I was asking recently about converting contour vectors to rasters and creating a DEM from the raster contours using r.surf.contour. My supplied vector contour layer has 20m intervals between contour lines. When I convert this to raster, the smallest region resolution I can use is 5m, otherwise v.to.rast runs out of memory. When I use r.surf.contour to create the DEM from the raster, which region resolution makes sense then? I set it to 20m because of the contour interval, but I'm starting to think that's not the best option because the interpolation between two contour lines should be as fine as possible. Any thoughts? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.to.rast gives empty raster layer
Hi Maris, Thanks, your suggestion worked. With the region's resolution set to 5m, I got some output. Regards Hanlie 2011/5/2, Maris Nartiss maris@gmail.com: g.region vect=contours_20m Also 20x20m cells wouldn't be too large?!? Maris. 2011/5/2, Hanlie Pretorius hanlie.pretor...@gmail.com: Hi, I'm working in GRASS 6.4.1 on Windows XP. I want to use r.surf.contour to convert contour lines into a DEM. However, I first need to convert my vector contour lines into raster contour lines. My region is set to the contours_20m vector layer, with a 20m resolution. v.to.rast --overwrite input=contours_20m output=contours_20m use=attr column=HEIGHT Loading data... Reading features... Writing raster map... Converted areas: 0 of 0 Converted points/lines: 1593 of 1593 v.to.rast complete. r.info for the result is: | Layer:contours_20m@CDNGI_DEM Date: Mon May 02 14:18:00 2011 | Mapset: CDNGI_DEM Login of Creator: hanlie | Location: SA_Lo_29E | DataBase: F:\grassdata | Title:Labels ( contours_20m@CDNGI_DEM ) | Timestamp: none | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1772 | Columns: 1750 | Total Cells: 3101000 |Projection: Transverse Mercator |N: -3121229.69982072S: -3156672.30728643 Res: 20.00147148 |E: -32474.90929218W: -67483.27873573 Res: 20.00478254 | Range of data:min = 1620 max = 2460 | Data Source: |Vector Map: contours_20m@CDNGI_DEM in mapset CDNGI_DEM |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=contours_20m@CDNGI_DEM layer=1 type=point,line,ar\ |ea output=contours_20m@CDNGI_DEM use=attr column=HEIGHT value\ |=1 rows=4096 As you can see from the number of categories (0), there is no information in this raster file. Can someone perhaps help me figure out why this happens? 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
[GRASS-user] v.to.rast gives empty raster layer
Hi, I'm working in GRASS 6.4.1 on Windows XP. I want to use r.surf.contour to convert contour lines into a DEM. However, I first need to convert my vector contour lines into raster contour lines. My region is set to the contours_20m vector layer, with a 20m resolution. v.to.rast --overwrite input=contours_20m output=contours_20m use=attr column=HEIGHT Loading data... Reading features... Writing raster map... Converted areas: 0 of 0 Converted points/lines: 1593 of 1593 v.to.rast complete. r.info for the result is: | Layer:contours_20m@CDNGI_DEM Date: Mon May 02 14:18:00 2011 | Mapset: CDNGI_DEM Login of Creator: hanlie | Location: SA_Lo_29E | DataBase: F:\grassdata | Title:Labels ( contours_20m@CDNGI_DEM ) | Timestamp: none | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1772 | Columns: 1750 | Total Cells: 3101000 |Projection: Transverse Mercator |N: -3121229.69982072S: -3156672.30728643 Res: 20.00147148 |E: -32474.90929218W: -67483.27873573 Res: 20.00478254 | Range of data:min = 1620 max = 2460 | Data Source: |Vector Map: contours_20m@CDNGI_DEM in mapset CDNGI_DEM |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=contours_20m@CDNGI_DEM layer=1 type=point,line,ar\ |ea output=contours_20m@CDNGI_DEM use=attr column=HEIGHT value\ |=1 rows=4096 As you can see from the number of categories (0), there is no information in this raster file. Can someone perhaps help me figure out why this happens? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] One chamfered corner with v.buffer on an extents rectangle
Hi, I'm working in GRASS 6.4.0 and using v.buffer applied to an extents rectangle. For some reason v.buffer is producing one chamfered corner in the buffered vector. The command I used was: v.buffer -s -c --overwrite input=c83a_extents@PERMANENT output=c83a_extents_5km_buffer@20m_DEM distance=5000. Buffering lines... Buffering areas... Building topology for vector map c83a_extents_5km_buffer... Registering primitives... 1 primitives registered 12 vertices registered Number of nodes: 1 Number of primitives: 1 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 0 Number of areas: - Number of isles: - Snapping boundaries... Breaking boundaries... Removing duplicates... Attaching islands... Building topology for vector map c83a_extents_5km_buffer... Building areas... 1 areas built 1 isles built Attaching islands... Number of nodes: 1 Number of primitives: 1 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 0 Number of areas: 1 Number of isles: 1 Number of areas without centroid: 1 Attaching centroids... Building topology for vector map c83a_extents_5km_buffer... Attaching centroids... Number of nodes: 2 Number of primitives: 2 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 1 Number of areas: 1 Number of isles: 1 Building topology for vector map c83a_extents_5km_buffer... Registering primitives... 2 primitives registered 13 vertices registered Building areas... 1 areas built 1 isles built Attaching islands... Attaching centroids... Number of nodes: 2 Number of primitives: 2 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 1 Number of areas: 1 Number of isles: 1 (Wed Apr 20 14:08:46 2011) Command finished (1 sec) I have placed a screenshot of the outcome on the web http://www.nedbib.za.net/v.buffer.png. The region was set to the DEM raster in the background. The black triangle is the extents and the red one the buffered polygon. Does someone know why v.buffer would create such a chamfered corner. I've played around with selecting different feature types and adding a minor distance value, but I always get the same result. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] One chamfered corner with v.buffer on an extents rectangle
Ok, thanks. The other three corners were rounded, but I eliminated the rounding with the -s option. Hanlie 2011/4/20, Markus Metz markus.metz.gisw...@googlemail.com: On Wed, Apr 20, 2011 at 2:21 PM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I'm working in GRASS 6.4.0 and using v.buffer applied to an extents rectangle. For some reason v.buffer is producing one chamfered corner in the buffered vector. The bug concerning the southwest corner is already fixed in 6.5, awaiting ok for backport to 6.4.2. The 3 other corners are also incorrect, I guess they should be rounded. At least they are with other areas of similar shape. Markus M The command I used was: v.buffer -s -c --overwrite input=c83a_extents@PERMANENT output=c83a_extents_5km_buffer@20m_DEM distance=5000. Buffering lines... Buffering areas... Building topology for vector map c83a_extents_5km_buffer... Registering primitives... 1 primitives registered 12 vertices registered Number of nodes: 1 Number of primitives: 1 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 0 Number of areas: - Number of isles: - Snapping boundaries... Breaking boundaries... Removing duplicates... Attaching islands... Building topology for vector map c83a_extents_5km_buffer... Building areas... 1 areas built 1 isles built Attaching islands... Number of nodes: 1 Number of primitives: 1 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 0 Number of areas: 1 Number of isles: 1 Number of areas without centroid: 1 Attaching centroids... Building topology for vector map c83a_extents_5km_buffer... Attaching centroids... Number of nodes: 2 Number of primitives: 2 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 1 Number of areas: 1 Number of isles: 1 Building topology for vector map c83a_extents_5km_buffer... Registering primitives... 2 primitives registered 13 vertices registered Building areas... 1 areas built 1 isles built Attaching islands... Attaching centroids... Number of nodes: 2 Number of primitives: 2 Number of points: 0 Number of lines: 0 Number of boundaries: 1 Number of centroids: 1 Number of areas: 1 Number of isles: 1 (Wed Apr 20 14:08:46 2011) Command finished (1 sec) I have placed a screenshot of the outcome on the web http://www.nedbib.za.net/v.buffer.png. The region was set to the DEM raster in the background. The black triangle is the extents and the red one the buffered polygon. Does someone know why v.buffer would create such a chamfered corner. I've played around with selecting different feature types and adding a minor distance value, but I always get the same result. 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
[GRASS-user] Installing r.inund.fluv on Windows
Hi, I'm working with GRASS 6.4.0 under QGIS in Windows XP. I'm posting here because I got no responses on the GRASS Windows email list. I'm trying to install the add-on r.inund.fluv (https://svn.osgeo.org/grass/grass-addons/raster/r.inund.fluv/). The readme file seems aimed at Linux users, so I followed the instructions that another user got from the author at http://gis.stackexchange.com/questions/6791/how-to-compile-in-fortran-in-order-to-use-an-add-on-in-grass. I managed to carry out all the steps in the above post and it created .exe files in the fortran_code folder. When I try to run r.inund.fluv or r.inund.fluv.bat from GRASS, I get the error r.inund.fluv is not recognized as an internal or external command. Perhaps someone can help me with the following questions: 1. What is the purpose of the .bat file? Does one actually run that? 2. What happened to the other steps in the make file? 3. Is there some way to run a makefile in Windows so that it will carry out all the instructions in that file? Any other ideas? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Assign attributes of start and enpoints to connecting lines
Hi Richard, Your suggestion worked in principle, but the whole operation needs more work. I ended up with the same node names in the start and end node columns, so I think a script is needed to identify the start and end nodes based on the direction of the flow of the river. An elevation column in the points layer could provide this information. Regards Hanlie Date: Sun, 12 Dec 2010 07:33:41 +1100 From: Richard Chirgwin rchirg...@ozemail.com.au Subject: Re: [GRASS-user] Assign attributes of start and enpoints to, connecting lines To: grass-user@lists.osgeo.org Message-ID: 4d03e025.7040...@ozemail.com.au Content-Type: text/plain; charset=ISO-8859-1; format=flowed Convert the start and end nodes to points using v.to.points, and then use v.distance? Richard ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Assign attributes of start and enpoints to connecting lines
Hi, I'm working in GRASS 6.4.0 on Windows XP. I'm creating a hydrological model and am using GIS to prepare the data. I've split my river layer into about 100 segments and I've digitised the start and end point of each river segment on a points layer. For the points layer, I calculate a 'name' from a prefix and a sequental number. What I would like to do, is to assign for each river segment the name of the start node and the name of the end node to the attribute table of the river segment layer. Can anyone give me pointers on how to do this? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Assign attributes of start and enpoints to connecting lines
Hi Micha, Thanks for the reply. I tried your method using the GUI, and it seems one cannot have a from_type of 'line', only of 'point' or 'centroid'. I tried: v.distance from=c83a_rivers_5km to=c83a_rivers_5km_junctions from_type=centroid to_type=point upload=to_attr column=Jname to_column=end_node And got: ERROR: Column Jname not found in table c83a_rivers_5km Jname is a column in the junctions layer. If I change the last two options around (obviously wrong): v.distance from=c83a_rivers_5km to=c83a_rivers_5km_junctions from_type=centroid to_type=point upload=to_attr column=end_node to_column=Jname I don't get an error but no updates either. I've found another (proprietary) program that can do this, but I was curious to see if it's possible in GRASS. Regards Hanlie 2010/12/9, Micha Silver mi...@arava.co.il: On 09/12/2010 16:38, Hanlie Pretorius wrote: Hi, I'm working in GRASS 6.4.0 on Windows XP. I'm creating a hydrological model and am using GIS to prepare the data. I've split my river layer into about 100 segments and I've digitised the start and end point of each river segment on a points layer. For the points layer, I calculate a 'name' from a prefix and a sequental number. What I would like to do, is to assign for each river segment the name of the start node and the name of the end node to the attribute table of the river segment layer. Can anyone give me pointers on how to do this? I think that v.distance can do this. begin by adding two attrib columns to your river layer, one for the start node label, and one for the end node label. Then run v.distance twice to get the labels from the 'name' column from each of the node layers. Finally concatenate the two labels together. So... # Add columns v.db.addcol rivers col=start_node varchar(8), end_node varchar(8), label varchar(16) # Use whatever size strings you need v.distance from=rivers to=end_nodes from_type=line to_type=point upload=to_attrib column=name to_col=end_node # and again for the start_nodes #Now merge the columns echo UPDATE rivers SET label=(start_node + end_node) | db.execute # I'm not sure the above '+' will work on all database connections. With dbf probably not... 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 http://www.surfaces.co.il/ Arava Development Co. +972-52-3665918 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] (no subject)
Hi, I'm running GRASS through the QGIS (v1.0.5) interface on Windows XP. The GRASS folder in the QGIS installation says 'grass-6.4.0svn'. I have a polygon layer with 22 polygons and a raster slope layer. When I run v.rast.stats with these two layers, two of the polygons have NULL values in the columns that v.rast.stats adds to the vector layer. These polygons are not the smallest or on the edges and I've used this procedure successfully before. Below are the region, v.info, r.info v.rast.stats outputs. Can anyone perhaps see why this happens? g.region -p projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -3122190 south: -3155760 west: -66510 east: -33390 nsres: 90 ewres: 90 rows: 373 cols: 368 cells: 137264 v.info map=subs22_clean...@chi ++ | Layer: subs22_clean...@chi | Mapset: CHI | Location:c83a_tm_29deg_e | Database:F:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: hanlie | Organization: | Source date: Wed Nov 17 13:24:57 2010 || | Type of Map: vector (level: 2) | | Number of points: 0 Number of areas: 22 | Number of lines:0 Number of islands:1 | Number of boundaries: 822 Number of faces: 0 | Number of centroids:22 Number of kernels:0 | | Map is 3D: No | Number of dblinks: 1 | | Projection: Transverse Mercator | N: -3122280S: -3155670 | E:-33570W:-66510 | | Digitization threshold: 0 | Comments: | ++ r.info map=c83a_dem_srtm_90m_sl...@permanent ++ | Layer:c83a_dem_srtm_90m_sl...@perma Date: Wed Nov 03 16:21:11 2010 | Mapset: PERMANENT Login of Creator: Administrator | Location: c83a_tm_29deg_e | DataBase: F:\Hanlie\grassdata | Title:percent slope ( c83a_dem_srtm_90m_slope ) | Timestamp: none |--- | | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 373 | Columns: 368 | Total Cells: 137264 |Projection: Transverse Mercator |N: -3122190S: -3155760 Res:90 |E: -33390W: -66510 Res:90 | Range of data:min = 0.023858 max = 110.669037 | | Data Source: |raster elevation file c83a_dem_srtm_...@permanent | | | Data Description: |generated by r.slope.aspect | | Comments: |slope map elev = c83a_dem_srtm_...@permanent |zfactor = 1.00 format = percent |min_slp_allowed = 0.00 | ++ v.rast.stats raster=c83a_dem_srtm_90m_sl...@permanent vector=subs22_clean...@chi =area =1 colprefix=slope Statistics calculated from raster map and uploaded to attribute table of vector map . Done. Successfully finished Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Euclidean distance between lines
Hi, I need to calculate the maximum Euclidean distance between lines in a vector layer. In ArcGIS I use Spatial Analyst to achieve this; it has a Euclidean distance tool that takes a line vector layer and produces a raster output, either as table or as an image. The output contains several statistics, including the maximum. I can't think of a way to do this in GRASS. v.dist gives only the minimum distance. r.distance also works with the shortest distance only. Can anyone perhaps suggest a method to accomplish this? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Conversion from shapefile to GRASS and back
Hi, I've imported a shapefile with 14 polygons using v.in.ogr: - v.in.ogr dsn=/media/0847147784/data/dwaf/liebenbergsvlei/myne/sa_tm_29_deg_east/C83A/c83a_14_subcatchments.shp output=test Datum unknown not recognised by GRASS and no parameters found Projection of input dataset and current location appear to match Layer: c83a_14_subcatchments Importing map 14 features... - Building topology for vector map test... Registering primitives... 14 primitives registered 16578 vertices registered Building areas... 7 areas built 4 isles built Attaching islands... Attaching centroids... Number of nodes: 9 Number of primitives: 14 Number of points: 0 Number of lines: 0 Number of boundaries: 14 Number of centroids: 0 Number of areas: 7 Number of isles: 4 Number of incorrect boundaries: 10 Number of areas without centroid: 7 - Cleaning polygons, result is not guaranteed! Building topology for vector map test... Number of nodes: 9 Number of primitives: 14 Number of points: 0 Number of lines: 0 Number of boundaries: 14 Number of centroids: 0 Number of areas: - Number of isles: - - Break polygons: - Remove duplicates: - Break boundaries: - Remove duplicates: - Clean boundaries at nodes: - Break boundaries: - Remove duplicates: - Clean boundaries at nodes: - Change dangles to lines: - Remove bridges: - Building topology for vector map test... Building areas... Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored 308 areas built 13 isles built Attaching islands... Number of nodes: 1221 Number of primitives: 3277 Number of points: 0 Number of lines: 0 Number of boundaries: 3277 Number of centroids: 0 Number of areas: 308 Number of isles: 13 Number of incorrect boundaries: 17 Number of areas without centroid: 308 Cannot calculate area centroid Cannot calculate area centroid Layer: c83a_14_subcatchments - - Building topology for vector map test... Registering primitives... 10001774 primitives registered 15588 vertices registered Building areas... Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored Area of size = 0.0 ignored 308 areas built 13 isles built Attaching islands... Attaching centroids... Number of nodes: 1447 Number of primitives: 1774 Number of points: 0 Number of lines: 0 Number of boundaries: 1539 Number of centroids: 235 Number of areas: 308 Number of isles: 13 Number of incorrect boundaries: 17 Number of duplicate centroids: 29 Number of areas without centroid: 102 - 196 areas represent more (overlapping) features, because polygons overlap in input layer(s). Such areas are linked to more than 1 row in attribute table. The number of features for those areas is stored as category in layer 2 14 input polygons Total area: 7.449545e+08 (308 areas) Overlapping area: 2.001273e-05 (196 areas) Area without category: 5.083098e-06 (71 areas) (Wed Nov 10 16:45:27 2010) Command finished (11 sec) - When I export the file to SHP again, it now contains 703 polygons: - v.out.ogr -e input=t...@permanent type=area dsn=/home/hanlie/ olayer=test The map contains islands. To preserve them in the output map, use the -c flag Exporting 308 areas (may take some
Re: [GRASS-user] Area weighting for vectors
2010/11/3, Micha Silver mi...@arava.co.il: On 11/03/2010 04:55 PM, Hanlie Pretorius wrote: Hi, I need to add some attributes from a soil polygon layer to a catchment polygon layer using area weighting. In other words, the contribution of the soil layer needs to be weighted by the amount with which its polygons overlap the catchment polygons. For example, suppose I want to weight soil infiltration rates. 30% of a catchment polygon overlaps with soil type 1 and 70% overlaps with soil type 2. Now I want a new attribute in the catchment polygon that takes 30% of the infiltration rate of soil type 1 and adds it to 70% of the infiltration rate of soil type 2. How can I do this in GRASS? The most straight forward way is using v.rast.stats. This module adds to a vector attribute table the univariate statistics from a raster. So you'll first convert the soil map to raster, then run v.rast.stats using the catchment polygons as the vector and the soil infiltration as the raster. Thanks, this worked. Is there a reason (other than no one has had time to do it) that a similar function doesn't exist for two vector layers? ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Area weighting for vectors
Hi, I need to add some attributes from a soil polygon layer to a catchment polygon layer using area weighting. In other words, the contribution of the soil layer needs to be weighted by the amount with which its polygons overlap the catchment polygons. For example, suppose I want to weight soil infiltration rates. 30% of a catchment polygon overlaps with soil type 1 and 70% overlaps with soil type 2. Now I want a new attribute in the catchment polygon that takes 30% of the infiltration rate of soil type 1 and adds it to 70% of the infiltration rate of soil type 2. How can I do this in GRASS? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.reclass with floating point result
2010/10/19, Glynn Clements gl...@gclements.plus.com: Hanlie Pretorius wrote: So, I have two problems here: 1) I can't get r.reclass to output an FCELL raster even though its help file suggest that only the input file needs to be a CELL raster. A reclass table maps integers to integers. You cannot generate a floating-point map by reclassing. Could you perhaps suggest a method to do so? I tried r.mapcalc, but I also got a conversion to integers. Somehow the default raster is an integer raster? -- 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] r.reclass with floating point result
Hi, I'm trying to reclassify a raster with these categories: |#|description |- |1|Sandy Loam |2|Sandy Clay Loam |*|no data To this: |#|description |- |1|32.7 |2|3 |*|no data So, I use a text file with the followin rules: 1 = 32.7 2 = 3 And the command r.reclass --overwrite input=c83a_soils output=c83a_soils_infiltration_mm_per_3_hours rules=infiltration_reclass_mm3hr.txt title=Saturated infiltration rate in mm per 3 hours 33.20 rounded up to 33 (Tue Oct 19 15:46:10 2010) Command finished (1 sec) The result is a raster with: | #|description |- | 3| |33| | *|no data I have also tried a rule file that says: 1 = 32.7 Sandy Loam 2 = 3 Sandy Clay Loam With the result | #|description |- | 3|Sandy Clay Loam |33|Sandy Loam | *|no data So, I have two problems here: 1) I can't get r.reclass to output an FCELL raster even though its help file suggest that only the input file needs to be a CELL raster. 2) I don't understand how r.reclass uses the rules file. Can someone perhaps help me? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[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
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
[GRASS-user] Script for r.proj and g.list issues
Hi, I'm working in GRASS 6.4.0RC6 and Ubuntu 10.04. I'm trying to write a Python script to project hundreds of files from one location to another. Since r.proj -l outputs text that will have to be parsed, I decided to rather follow the method explained at this webpage (http://www.mail-archive.com/grass-...@lists.osgeo.org/msg14180.html), namely to export a list of rasters to a text file using g.list. If I run g.list without the mapset option, I get a list of all the rasters, but I'm particularly interested in the mapset called '2000_02_february_trmm'. I include some of the output below to prove that the mapset contains some rasters: - g.list type=rast -- raster files available in mapset 2000_02_february_trmm: 2201_00_3B42RT 2208_09_3B42RT 2215_18_3B42RT ... 2208_06_3B42RT 2215_15_3B42RT 2223_00_3B42RT raster files available in mapset PERMANENT: MCD12Q1.A2001001.005_land_cover_type1 ... test_border_nulls -- (Mon Sep 13 12:22:14 2010) Command finished (0 sec) - However, g.list is behaving strangely when I ask it to list the rasters in mapset 2000_02_february_trmm: - g.list type=rast mapset=2000_02_february_trmm -- no raster files available in mapset 2000_02_february_trmm -- (Mon Sep 13 12:23:11 2010) Command finished (0 sec) -- So, I have two questions: 1) Am I following the right method to project the rasters from one location to the other, or does an easier way exist? 2) Why is g.list not finding my mapset when it clearly exists? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Script for r.proj and g.list issues
2010/9/13, Sylvain Maillard sylvain.maill...@gmail.com: rasters in mapset 2000_02_february_trmm: - g.list type=rast mapset=2000_02_february_trmm -- no raster files available in mapset 2000_02_february_trmm -- (Mon Sep 13 12:23:11 2010) Command finished (0 sec) -- you made a mistake in the name of the mapset: you used 2000_02_february_trmm, try again without a at the end ... Sylvain My bad, sorry! ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Adding options to grass.run_command
Hi, I'm writing my first GRASS Python script and am having trouble adding options to commands that I want to run. For example, the following works: - grass.run_command(r.in.xyz, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) - but this doesn't: - grass.run_command(r.in.xyz, --overwrite, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) - neither does this: - grass.run_command(r.in.xyz --overwrite, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) - How can I get this to work, or is there a better way of doing it using a command other than grass.run_command? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Adding options to grass.run_command
Thanks :-) 2010/9/3, Martin Landa landa.mar...@gmail.com: 2010/9/3 Hanlie Pretorius hanlie.pretor...@gmail.com: grass.run_command(r.in.xyz, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) better `input = os.path.join(directory, file)` but this doesn't: - grass.run_command(r.in.xyz, --overwrite, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) neither does this: - grass.run_command(r.in.xyz --overwrite, input=directory+'/'+file, output=rainraster, fs=',', x=2, y=1) should be `overwrite = True` instead `--overwrite` Martin -- Martin Landa landa.martin gmail.com * http://gama.fsv.cvut.cz/~landa ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Raster resolution problems
Hi, I'm using GRASS6.4SVN on Windows XP. I want to import TRMM precipitation data into a GRASS projected (Transverse Mercator) location. I have previously successfully import this data to a GCS (WGS84) location using the procedure at http://grass.osgeo.org/wiki/Import_XYZ. My problem is the grid resolution. The TRMM grid consists of 0.25° squares, but projected they have different sizes. Here is the output from v.report (area in km2) for a vector layer with the TRMM grid squares, projected into the Transverse Mercator location: - cat|LabelX|LabelY|LblOffsetX|LblOffsetY|Label|gRow|gColumn|RowCol|area 5|28.25|-28.75|-20|-20|28.25677.428350142518 11|28.5|-28.5|-20|-20|28.5678.961903132767 6|28.25|-28.5|-20|-20|28.25679.002141385312 1|28|-28.5|-20|-20|28679.062503849138 12|28.5|-28.25|-20|-20|28.5680.522519692921 7|28.25|-28.25|-20|-20|28.25680.563040996854 2|28|-28.25|-20|-20|28680.623828116383 13|28.5|-28|-20|-20|28.5682.070220966439 8|28.25|-28|-20|-20|28.25682.111024322052 3|28|-28|-20|-20|28682.172234597964 14|28.5|-27.75|-20|-20|28.5683.604982570918 9|28.25|-27.75|-20|-20|28.25683.646066944101 4|28|-27.75|-20|-20|28683.70769882544 15|28.5|-27.5|-20|-20|28.5685.126780361735 10|28.25|-27.5|-20|-20|28.25685.168144684034 - The projected grid 'squares' vary in size from 677km2 to 685km2. So how am I supposed to set the resolution of the region when I do the import? (I previously tried to reproject the TRMM data from the GCS location to the Transverse Mercator location, but then I ran into the same problem - resolution. Any ideas? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Compiling GRASS with debugging symbols
Hi I would like to compile GRASS with extra debugging information and apparently one needs to enable 'debugging symbols'. Which of the configuration options would add debugging symbols to the compilation? I can't see any in the list produced by - ./configure --help - that obviously include this option. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Python scripts: 'command not found'
2010/8/19, razmjoo...@faunalia.co.uk razmjoo...@faunalia.co.uk: In Ubuntu, you can right click on the file and tick the box to make it executable Right Click Properties Permissions Allow executing files as program Cheers Sab Thanks, tried this and it workedfor the script in the /usr/lib/grass64/scripts folder. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Python scripts: 'command not found'
Hi, I'm working in GRASS 6.4RC6 on Ubuntu 10.04 and I'm trying to learn how to script GRASS using Python. I have taken the example at http://grass.osgeo.org/programming6/pythonlib.html that checks if a vector is 3D and saved it to my home folder. When I try to run it though, GRASS can't find the file: - GRASS 6.4.0RC6 (sa_wgs84):~ ls *.py check_3d.py proxy.py GRASS 6.4.0RC6 (sa_wgs84):~ check_3d.py check_3d.py: command not found I also tried placing the script in the GRASS scripts folder using the name 'v.check3d' and got the following errors: - GRASS 6.4.0RC6 (sa_wgs84):~ v.check3d bash: /usr/lib/grass64/scripts/v.check3d: Permission denied GRASS 6.4.0RC6 (sa_wgs84):~ sudo v.check3d [sudo] password for hanlie: sudo: v.check3d: command not found - Can someone perhaps help me? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Python scripts: 'command not found'
2010/8/19, Nikos Alexandris nikos.alexand...@felis.uni-freiburg.de: Better not to use sudo (as a general advice). The script needs to be made executable for your user-name, i.e.: sudo chown hanlie:hanlie v.check3d Ran 'sudo chown hanlie:hanlie v.check3d' from the /usr/lib/grass64/scripts directory with no error messages, but GRASS still gives a permission error: - GRASS 6.4.0RC6 (sa_wgs84):~ v.check3d bash: /usr/lib/grass64/scripts/v.check3d: Permission denied - ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.sim.water buffer overflow
Hi, I'm using GRASS64 RC6 on Linux 10.04 (64-bit) and I'm trying to run r.sim.water on a catchment about 5000 km2 in size. I started with a very fine resolution in my region (00:00:00.18) and got the following: - GRASS 6.4.0RC6 (world_wgs84):~ r.sim.water elevin=lieb_dem_25m_clipped dxin=lieb_dem_25m_clipped_dx dyin=lieb_dem_25m_clipped_dy rain_val=2.0 manin_val=0.1 infil_val=13.0 depth=rsimwater.depth.18s disch=rsimwater.disch.18s default nwalk=600505650, rwalk=600505650.00 Killed GRASS 6.4.0RC6 (world_wgs84):~ (python:1572): Pango-CRITICAL **: pango_layout_get_cursor_pos: assertion `index = 0 index = layout-length' failed ^C Since this resolution gives about 300 million cells in the region, I progressively coarsened the resolution, but I get buffer overflows even with just 2736 cells in the region (resolution = 00:01:00) Below is an example of the error I got at a resolution of 1 minute: - GRASS 6.4.0RC6 (world_wgs84):~ r.sim.water elevin=lieb_dem_25m_clipped dxin=lieb_dem_25m_clipped_dx dyin=lieb_dem_25m_clipped_dy rain_val=50.0 manin_val=0.1 infil_val=13.0 depth=rsimwater.depth_1min disch=rsimwater.disch_1min --overwrite default nwalk=5472, rwalk=5472.00 Min elevation = 1500.00 m Max elevation = 2351.17 m Mean Source Rate (rainf. excess or sediment) = 0.06 m/s or kg/m2s Mean flow velocity = 1.135392 m/s Mean Mannings = 0.170467 Number of iterations = 41128 cells Time step = 0.00 s *** buffer overflow detected ***: r.sim.water terminated === Backtrace: = /lib/libc.so.6(__fortify_fail+0x37)[0x7ffab34e7207] /lib/libc.so.6(+0xfe0c0)[0x7ffab34e60c0] /lib/libc.so.6(+0xfd529)[0x7ffab34e5529] /lib/libc.so.6(_IO_default_xsputn+0xcc)[0x7ffab345dd1c] /lib/libc.so.6(_IO_vfprintf+0x3d34)[0x7ffab34310d4] /lib/libc.so.6(__vsprintf_chk+0x99)[0x7ffab34e55c9] /lib/libc.so.6(__sprintf_chk+0x7f)[0x7ffab34e550f] /usr/lib/grass64/lib/libgrass_sim.so(output_data+0x1b1a)[0x7ffab5f32d35] r.sim.water(main+0x1041)[0x4037e5] /lib/libc.so.6(__libc_start_main+0xfd)[0x7ffab3406c4d] r.sim.water[0x4026e9] === Memory map: 0040-00405000 r-xp 08:05 545723 /usr/lib/grass64/bin/r.sim.water 00604000-00605000 r--p 4000 08:05 545723 /usr/lib/grass64/bin/r.sim.water 00605000-00606000 rw-p 5000 08:05 545723 /usr/lib/grass64/bin/r.sim.water 00606000-12bc4000 rw-p 00:00 0 14a19000-14ab6000 rw-p 00:00 0 [heap] 7ffaaa3c1000-7ffaaa3cd000 r-xp 08:05 262845 /lib/libnss_files-2.11.1.so 7ffaaa3cd000-7ffaaa5cc000 ---p c000 08:05 262845 /lib/libnss_files-2.11.1.so 7ffaaa5cc000-7ffaaa5cd000 r--p b000 08:05 262845 /lib/libnss_files-2.11.1.so 7ffaaa5cd000-7ffaaa5ce000 rw-p c000 08:05 262845 /lib/libnss_files-2.11.1.so 7ffaaa5ce000-7ffaaa5d8000 r-xp 08:05 263050 /lib/libnss_nis-2.11.1.so 7ffaaa5d8000-7ffaaa7d7000 ---p a000 08:05 263050 /lib/libnss_nis-2.11.1.so 7ffaaa7d7000-7ffaaa7d8000 r--p 9000 08:05 263050 /lib/libnss_nis-2.11.1.so 7ffaaa7d8000-7ffaaa7d9000 rw-p a000 08:05 263050 /lib/libnss_nis-2.11.1.so 7ffaaa7d9000-7ffaaa7e1000 r-xp 08:05 262840 /lib/libnss_compat-2.11.1.so 7ffaaa7e1000-7ffaaa9e ---p 8000 08:05 262840 /lib/libnss_compat-2.11.1.so 7ffaaa9e-7ffaaa9e1000 r--p 7000 08:05 262840 /lib/libnss_compat-2.11.1.so 7ffaaa9e1000-7ffaaa9e2000 rw-p 8000 08:05 262840 /lib/libnss_compat-2.11.1.so 7ffaaa9e2000-7ffaaa9e5000 r-xp 08:05 261719 /lib/libgpg-error.so.0.4.0 7ffaaa9e5000-7ffaaabe4000 ---p 3000 08:05 261719 /lib/libgpg-error.so.0.4.0 7ffaaabe4000-7ffaaabe5000 r--p 2000 08:05 261719 /lib/libgpg-error.so.0.4.0 7ffaaabe5000-7ffaaabe6000 rw-p 3000 08:05 261719 /lib/libgpg-error.so.0.4.0 7ffaaabe6000-7ffaaabf6000 r-xp 08:05 137210 /usr/lib/libtasn1.so.3.1.7 7ffaaabf6000-7ffaaadf5000 ---p 0001 08:05 137210 /usr/lib/libtasn1.so.3.1.7 7ffaaadf5000-7ffaaadf6000 r--p f000 08:05 137210 /usr/lib/libtasn1.so.3.1.7 7ffaaadf6000-7ffaaadf7000 rw-p 0001 08:05 137210 /usr/lib/libtasn1.so.3.1.7 7ffaaadf7000-7ffaaae1 r-xp 08:05 137156 /usr/lib/libsasl2.so.2.0.23 7ffaaae1-7ffaab00f000 ---p 00019000 08:05 137156 /usr/lib/libsasl2.so.2.0.23 7ffaab00f000-7ffaab01 r--p 00018000 08:05 137156 /usr/lib/libsasl2.so.2.0.23 7ffaab01-7ffaab011000 rw-p 00019000 08:05 137156 /usr/lib/libsasl2.so.2.0.23 7ffaab011000-7ffaab027000 r-xp 08:05 263380 /lib/libresolv-2.11.1.so 7ffaab027000-7ffaab226000 ---p 00016000 08:05 263380 /lib/libresolv-2.11.1.so 7ffaab226000-7ffaab227000 r--p 00015000 08:05 263380 /lib/libresolv-2.11.1.so 7ffaab227000-7ffaab228000 rw-p 00016000 08:05 263380 /lib/libresolv-2.11.1.so 7ffaab228000-7ffaab22a000 rw-p 00:00 0 7ffaab22a000-7ffaab22c000 r-xp 08:05 261726 /lib/libkeyutils-1.2.so 7ffaab22c000-7ffaab42b000 ---p 2000 08:05 261726 /lib/libkeyutils-1.2.so 7ffaab42b000-7ffaab42c000 r--p 1000
[GRASS-user] Re: Importing multiple files with r.in.xyz
Thanks for the reply. I'll try that if I don't manage with the Linux version. However, previously when I worked with the SVN version, I was advised to rather use the stable version. 2010/8/10, Helmut Kudrnovsky hel...@web.de: Hi, I'm working with WinGRASS 6.4RC6 on Win XP [...] please try also the latest nightly builds of WinGrass64svn (http://josef.fsv.cvut.cz/wingrass/grass64/), because release candidate 6 of Grass64 is now a little bit outdated. best regards Helmut ___ Neu: WEB.DE De-Mail - Einfach wie E-Mail, sicher wie ein Brief! Jetzt De-Mail-Adresse reservieren: https://produkte.web.de/go/demail02 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Re: Importing multiple files with r.in.xyz
Had a look at the article. Interesting process. However, the TRMM rainfall data is already spatially distributed, albeit at a coarse resolution (0.25 degrees). 2010/8/11, razmjoo...@faunalia.co.uk razmjoo...@faunalia.co.uk: Hanlie You might find this useful: http://www.surfaces.co.il/?p=578 It's in qgis, but still you can use GRASS plugin to interact with your datasets. Cheers Saber Thanks for the reply. I'll try that if I don't manage with the Linux version. However, previously when I worked with the SVN version, I was advised to rather use the stable version. 2010/8/10, Helmut Kudrnovsky hel...@web.de: Hi, I'm working with WinGRASS 6.4RC6 on Win XP [...] please try also the latest nightly builds of WinGrass64svn (http://josef.fsv.cvut.cz/wingrass/grass64/), because release candidate 6 of Grass64 is now a little bit outdated. best regards Helmut ___ Neu: WEB.DE De-Mail - Einfach wie E-Mail, sicher wie ein Brief! Jetzt De-Mail-Adresse reservieren: https://produkte.web.de/go/demail02 ___ 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] Importing multiple files with r.in.xyz
My apologies for my ignorance, but how do I get this to run in Linux? Am using Ubuntu 10.04 and GRASS 6.4RC6. I entered the following into the GRASS text window and am testing it with a list of two files (I removed the region settings because the files all have the same region, which I have already set): - GRASS 6.4.0RC6 (world_wgs84):~ cat lieb_files_test.txt 3B42.000201.12.6.nc.lieb.txt 3B42.000201.15.6.nc.lieb.txt GRASS 6.4.0RC6 (world_wgs84):~ cat lieb_files_test.txt | while read line; do echo $line+ being processed raster=${$line:(-12)} # cut .nc.lieb.txt from raster name r.in.xyz --overwrite input=$line output=raster method=mean type=FCELL fs=, x=2 y=1 z=3 - As you can see, when I press enter after the r.in.xyz line, I just get another prompt. Am I supposed to save it to a separate file and run this as a script? If so, where should I put this script? 2010/8/10, Saber Razmjooei razmjoo...@faunalia.co.uk: cat list_of_files.txt | while read line; do echo $line+Being processed # or whaterver you want to do with the $line variable raster=$line r.in.xyz -s -g input=$line output=$line tmpRegion myregion= `head -n 1 tmpRegion` g.region $myregion r.in.xyz --overwrite input=$line output=$line fs=, done and your list_of_files.txt is raster1 raster2 . Hope that helps Saber ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing multiple files with r.in.xyz
2010/8/11, Vincent Bain b...@toraval.fr: Hi, short intrusion in this current thread, without knowing exactly what it is related to : I only mean to suggest you to close the while loop... This prompt means the shell expects the end of your statement ; here you should type done to close the while instruction. Adding 'done' worked and I got the script running from the command line. Thanks to everyone who chipped in. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing multiple files with r.in.xyz
2010/8/11, Micha Silver mi...@arava.co.il: Samber's method will surely work, but you might more simply try as follows: # start grass in a location which matches the txt file data # change to the directory where your txt files are, then do for infile in *.txt; do outrast=`basename ${infile} .txt`; r.in.xyz in=${infile} out=${outrast} fs=,; done This works nicely, except that my input files have names like: 3B42.000201.0.6.nc.lieb.txt 3B42.000201.3.6.nc.lieb.txt and I want the output rasters to be named 3B42.000201.0.6 3B42.000201.3.6 In fact, it would be best if they could be named: 0201.0.6 0201.3.6 I'm struggling with this string manipulation in the shell scripting language. # that's it ;-) -- Micha ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing multiple files with r.in.xyz
2010/8/11, Nikos Alexandris nikos.alexand...@felis.uni-freiburg.de: Micha Silver mi...@arava.co.il: Samber's method will surely work, but you might more simply try as follows: # start grass in a location which matches the txt file data # change to the directory where your txt files are, then do for infile in *.txt; do outrast=`basename ${infile} .txt`; r.in.xyz in=${infile} out=${outrast} fs=,; done Hanlie Pretorius wrote: This works nicely, except that my input files have names like: 3B42.000201.0.6.nc.lieb.txt 3B42.000201.3.6.nc.lieb.txt and I want the output rasters to be named 3B42.000201.0.6 3B42.000201.3.6 In fact, it would be best if they could be named: 0201.0.6 0201.3.6 I'm struggling with this string manipulation in the shell scripting language. Try the following (note the backticks in the beginning before echo and in the end as well): output=`echo ${infile} | cut -d. -f2,3,4 | sed 's/^00//'` Yes, this worked :-) Thanks. Nikos ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Importing multiple files with r.in.xyz
Hi, I'm working with WinGRASS 6.4RC6 on Win XP. I'm following the steps at http://grass.osgeo.org/wiki/Import_XYZ to import hundreds files containing TRMM rainfall data. When I try to follow the instruction: - Just amend above procedure to use wildcards. Change in above example all occurencies of cat VTL2733.XYZ | ... to cat *.XYZ | ... and use a more reasonable output name of course. That's all to import even thousands of files (tiled DEM) easily. - Not sure what a 'reasonable output name is', so I run it as shown below and get the following error: - cat F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\*.txt| r.in.xyz -s in=- fs=, out=test cat: F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\*.txt: No such file or directory - If I try it with one file only, it works: - cat F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\3B42.000201.0.6.nc.lieb.txt| r.in.xyz -s in=- fs=, out=test Range: x:-28.625000-27.375000 y: 28.125000 28.625000 z: 0.00 0.00 - How do I get this to work? And how would I get the output filenames to be related to the input filenames? For example, if the input filename is 3B42.000201.0.6.nc.lieb.txt, then I want the output raster to be named '3B42.000201.0.6'. Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Clipping rasters
Hi, I thought I created a clipped DEM by rasterising an irregular vector boundary and then using a mask to create the clipped DEM from a rectangular DEM. For example: - # convert vector border to raster v.to.rast input=border output=border use=cat type=area layer=1 value=1 rows=4096| # set the mask to the raster border r.mask border # create the clipped DEM r.mapcalc dem_clipped=dem_full # remove the mask r.mask -r - But it turns out I was wrong, the 'clipped' raster just contains NULL values outside of the border raster. What did I do wrong? Is there a better way to clip a raster to a vector boundary? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing multiple files with r.in.xyz
Thanks for the suggestions. I think I'm going to try the Linux version as soon as I can get to a Linux machine. 2010/8/10, Saber Razmjooei razmjoo...@faunalia.co.uk: Hanile Try this: 1- Install Paths: http://www.textpad.com/add-ons/files/utilities/paths.zip It's a useful tool to copy the path and filenames in windows. It adds a function to rightclick: Pathcopy. You can select all your rasters and rightclick to copy the paths: Paste them in a text file c:\example\raster1.xyz c:\example\raster2.xyz ... 2- change them to something like this (use find/replace..macros, etc:) r.in.xyz -s -g input=c:\example\raster1.xyz output=raster1 tmpRegion SET /p myregion= tmpRegion DEL tmpRegion g.region %myregion% r.in.xyz --overwrite input=c:\example\raster1.xyz output=raster1 fs=, 3- Save that as a import.bat 4- From the Start menu, run Grass in Txt mode, change directory to the folder where import.bat is 5- import.bat If you google around you should be able to make the script a bit smarter using DOS and loop (http://www.robvanderwoude.com/ntfor.php). In linux for example you can shrink it to something like this: cat list_of_files.txt | while read line; do echo $line+Being processed # or whaterver you want to do with the $line variable raster=$line r.in.xyz -s -g input=$line output=$line tmpRegion myregion= `head -n 1 tmpRegion` g.region $myregion r.in.xyz --overwrite input=$line output=$line fs=, done and your list_of_files.txt is raster1 raster2 . Hope that helps Saber On Tue, 2010-08-10 at 14:59 +0200, Hanlie Pretorius wrote: Hi, I'm working with WinGRASS 6.4RC6 on Win XP. I'm following the steps at http://grass.osgeo.org/wiki/Import_XYZ to import hundreds files containing TRMM rainfall data. When I try to follow the instruction: - Just amend above procedure to use wildcards. Change in above example all occurencies of cat VTL2733.XYZ | ... to cat *.XYZ | ... and use a more reasonable output name of course. That's all to import even thousands of files (tiled DEM) easily. - Not sure what a 'reasonable output name is', so I run it as shown below and get the following error: - cat F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\*.txt| r.in.xyz -s in=- fs=, out=test cat: F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\*.txt: No such file or directory - If I try it with one file only, it works: - cat F:\Hanlie\UCT\M.Sc.\Data\TRMM\2000\02_Februrarie\3B42.000201.0.6.nc.lieb.txt| r.in.xyz -s in=- fs=, out=test Range: x:-28.625000-27.375000 y: 28.125000 28.625000 z: 0.00 0.00 - How do I get this to work? And how would I get the output filenames to be related to the input filenames? For example, if the input filename is 3B42.000201.0.6.nc.lieb.txt, then I want the output raster to be named '3B42.000201.0.6'. 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] Clipping rasters
2010/8/10, Nikos Alexandris nikos.alexand...@felis.uni-freiburg.de: Hanlie Pretorius wrote: I thought I created a clipped DEM by rasterising an irregular vector boundary and then using a mask to create the clipped DEM from a rectangular DEM. For example: - #first, check and set your region of interest g.region vect=border -pa Tried this with the same result, as Jarek points out in the second reply. So, I assume it's impossible to have a raster with an irregular boundary in GRASS? # convert vector border to raster v.to.rast input=border output=border use=cat type=area layer=1 value=1 rows=4096| # set the mask to the raster border r.mask border # create the clipped DEM r.mapcalc dem_clipped=dem_full # remove the mask r.mask -r - But it turns out I was wrong, the 'clipped' raster just contains NULL values outside of the border raster. What did I do wrong? Is there a better way to clip a raster to a vector boundary? So, always check the region ;-). I think this is the best (and only?) way to clip. Regards, Nikos ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Filling Null values in a MODIS image
Hi, I have a MODIS land cover image that contains NULL values according to r.univar: - r.univar map=modis_land_cover_type1_2001_prim...@permanent total null and non-null cells: 41884 total null cells: 17145 - So, I tried to fill it with r.fillnulls and I get the following: - r.fillnulls input=mcd12q1.a2001001.005_land_cover_ty...@permanent output=MCD12Q1.A2001001.005_land_cover_type1_filled Locating and isolating NULL areas... Reading input raster map r_fillnulls_1...@permanent... Writing output raster map r_fillnulls_1911.buf... Creating interpolation points... Extracting points... Building topology for vector map vecttmp_fillnulls_1911... Registering primitives... 0 primitives registered 0 vertices registered Building areas... 0 areas built 0 isles built Attaching islands... Attaching centroids... Number of nodes: 0 Number of primitives: 0 Number of points: 0 Number of lines: 0 Number of boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of isles: 0 r.to.vect complete. Interpolating 0 points Not sufficient points to interpolate. Maybe no hole(s) to fill in the current map region? Removing raster MASK Removing raster r_fillnulls_1911 Removing raster r_fillnulls_1911.buf Removing raster r_fillnulls_1911_filled Raster map r_fillnulls_1911_filled not found r_fillnulls_1911_filled nothing removed Removing vector vecttmp_fillnulls_1911 (Wed Aug 4 10:18:22 2010) Command finished (0 sec) - Output of g.region -p: - g.region -p rast=modis_land_cover_type1_2001_prim...@permanent res=0.00423094 projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 27:19:08.833466S south: 28:30:53.475952S west: 28:03:57.840271E east: 28:41:37.727966E nsres: 0:00:15.210751 ewres: 0:00:15.269511 rows: 283 cols: 148 cells: 41884 - Output of r.info: - r.info map=modis_land_cover_type1_2001_prim...@permanent ++ | Layer:modis_land_cover_type1_2001_p Date: Thu Jul 29 15:17:55 2010 | Mapset: PERMANENT Login of Creator: hanlie | Location: world_wgs84 | DataBase: /media/0847147784/grassdata | Title: ( MCD12Q1.A2001001.005_land_cover_type1 ) | Timestamp: none | | | Type of Map: raster Number of Categories: 14 | Data Type:CELL | Rows: 11957 | Columns: 6277 | Total Cells: 75054089 |Projection: Latitude-Longitude |N: 27:19:08.833466SS: 28:30:53.475952S Res: 0:00:00.36001 |E: 28:41:37.727966EW: 28:03:57.840271E Res: 0:00:00.36002 | Range of data:min = 1 max = 14 | | Data Description: |generated by r.in.gdal | | Comments: |r.in.gdal input=F:\Hanlie\UCT\M.Sc\Data\modis\MOD12Q1\MCD12Q1.A2001\ |001.h20v11.005.2009342153812_land_cover_type1.tif output=MCD12Q1.A\ |2001001.005_land_cover_type1 | ++ I have tried to see the NULLS by displaying everything in white, except for NULLS displayed in red, but I didn't see anything. Does anyone know how I can find out if there are actually NULL value and where they are? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: Filling Null values in a MODIS image
2010/8/4, Hamish hamis...@yahoo.com: Hanlie wrote: | Rows: 11957 | Columns: 6277 ... I have tried to see the NULLS by displaying everything in white, except for NULLS displayed in red, but I didn't see anything. Does anyone know how I can find out if there are actually NULL value and where they are? try zooming in. Unless you have a 12000x7000 monitor, some cells will be hidden. extract with r.mapcalc and buffer with r.buffer if they are still hard to find. Thanks for all the replies. I first tried finding the NULLS using r.mapcalc, and it turns out that the NULLS are the cells in the GRASS region (rectangular), but outside my study region (irregular geographical border). I was under the misconception that g.region set to a raster file functions like a mask. About the r.neighbors procedure, why is it necessary to follow with r.patch? Surely the result of r.neighbors will be without NULLS? I did find, though, that r.neighbours resulted in raster with extended borders compared to the original, which makes sense. One would then have to clip it again to the study area. Or use a mask of the original area before running r.patch? Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Location of r.sim.water code
Hi, I work with GRASS 6.4RC6 and would like to look at the code for the r.sim.water module. I have been looking at http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_6_4/raster/, but I can't find it there. Can someone please tell me where I the code is kept? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Issues when reprojecting from a GCS to a PCS
Hi Everyone, I'm working in GRASS 6.4RC6 on Win XP and am trying to reproject TRMM rainfall data from its native GCS (WGS84) to a local PCS. I'm using the PCS because my other important data, specifically DEMs, are provided in it. TRMM data come in 0.25°x0.25° squares. My workfow is as follows: 1. Import the TRMM data into a WGS84 location using r.in.xyz. This gives a raster with the following information: ++ | Layer:3b42.00021...@permanentDate: Thu Jul 15 14:42:27 2010 | Mapset: PERMANENT Login of Creator: hanlie | Location: world_wgs84 | DataBase: F:\Hanlie\grassdata | Title:Raw x,y,z data binned into a raster grid by cell mean ( 3B42.000 | Timestamp: none |--- | | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 6 | Columns: 3 | Total Cells: 18 |Projection: Latitude-Longitude |N: 27:15SS: 28:45S Res: 0:15 |E: 28:45EW:28E Res: 0:15 | Range of data:min = 0.00 max = 0.070614 | | Data Source: |F:\Hanlie\UCT\M.Sc\Data\TRMM\2000\02_Februarie\3B42.000215.0.6.nc.lieb. | | | Data Description: |generated by r.in.xyz | | Comments: |r.in.xyz input=F:\Hanlie\UCT\M.Sc\Data\TRMM\2000\02_Februarie\3B42.\ |000215.0.6.nc.lieb.txt output=3B42.000215.0 method=mean type=F\ |CELL x=2 y=1 z=3 zscale=1.0 percent=100 | ++ 2. Set the region to this raster file and create a vector map from it using v.in.region. This gives a vector with the following information: ++ | Layer: lieb_trmm_reg...@permanent | Mapset: PERMANENT | Location:world_wgs84 | Database:F:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: hanlie | Organization: | Source date: Thu Jul 15 14:55:13 2010 || | Type of Map: vector (level: 2) | | Number of points: 0 Number of areas: 1 | Number of lines:0 Number of islands:1 | Number of boundaries: 1 Number of faces: 0 | Number of centroids:1 Number of kernels:0 | | Map is 3D: No | Number of dblinks: 0 | | Projection: Lat/Lon | N:27:15SS:28:45S | E:28:45EW: 28E | | Digitization threshold: 0 | Comments: | ++ 3. Change to the PCS location and reproject the TRMM region using v.proj. This gives a vector with the following information: ++ | Layer: lieb_trmm_reg...@permanent | Mapset: PERMANENT | Location:sa_tm_29deg_E | Database:F:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: hanlie | Organization: | Source date: Thu Jul 15 14:55:13 2010 || | Type of Map: vector (level: 2) | | Number of points: 0 Number of areas: 1 | Number of lines:0 Number of islands:1 | Number of boundaries: 1 Number of faces: 0 | Number of centroids:1 Number of kernels:0 | | Map is 3D: No | Number of dblinks: 0 | | Projection: Transverse Mercator | N: -3015356.39959204S: -3181970.90748457 | E:-24418.1519951W: -99037.39660484 | | Digitization threshold: 0 | Comments: | ++ 4. Set the region of the PCS to the TRMM region: g.region -p vect=lieb_trmm_reg...@permanent nsres=25 ewres=25 which gives: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -3015356.39959204 south: -3181970.90748457 west: -99037.39660484 east: -24418.1519951 nsres: 24.99842579 ewres: 24.9980719 rows: 6665 cols: 2985 cells: 19895025 I use this resolution because I don't really know what else to use. When I reproject the TRMM GRID into the PCS, I get polygons with differing areas, so I can't use the TRMM grid size as resolution: v.report map=lieb_trmm_squa...@permanent option=area units=kilometers Displaying column types/names for
Re: [GRASS-user] Re: Importing Dems with r.in.xyz
2010/6/5, Hamish hamis...@yahoo.com: Hanlie wrote: At this point, g.region reports 1146474 cells in the region, while I have 1146370 lines of coordinates in my file. ... So it looks like there are about 100 coordinates missing from the ASCII ASCII file. 0.01% .. Maybe holes in the data? perhaps this: https://trac.osgeo.org/grass/ticket/123 ?? I don't think it's this bug because this bug discards only one line of data. I don't get any data in because the number of coordinate pairs in the file is less than the number of cells in the defined region. I was thinking perhaps importing the points as vectors, converting them to raster and then doing a nearest neighbour or IDW interpolation to fill the gaps. At least then I'll be able to see where the gaps are and limit the interpolated pixels using a mask? No need to do anything different to find the missing pixels. Inspecting the output of r.univar with r.in.xyz's method=n maps can be very useful for troubleshooting. from the help page: Gridded data If data is known to be on a regular grid r.in.xyz can reconstruct the map perfectly as long as some care is taken to set up the region correctly and that the data's native map projection is used. A typical method would involve determining the grid resolution either by examining the data's associated documentation or by studying the text file. Next scan the data with r.in.xyz's -s (or -g) flag to find the input data's bounds. GRASS uses the cell-center raster convention where data points fall within the center of a cell, as opposed to the grid-node convention. Therefore you will need to grow the region out by half a cell in all directions beyond what the scan found in the file. After the region bounds and resolution are set cor- rectly with g.region, run r.in.xyz using the n method and verify that n=1 at all places. r.univar can help. Once you are confident that the region exactly matches the data proceed to run r.in.xyz using one of the mean, min, max, or median methods. With n=1 throughout, the result should be identical regardless of which of those methods are used. with the n map you might use r.mapcalc to extract the NULL cells as some value, then r.out.xyz or r.to.vect on th extracts to highlight where they are. Or maybe you get lucky with r.colors with nv set to bright magenta on the original data. Thanks, I'll try this to find where the holes in the data are. Hamish ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Re: Importing Dems with r.in.xyz
2010/6/7, Hamish hamis...@yahoo.com: Hanlie Pretorius wrote: I don't think it's this bug because this bug discards only one line of data. I don't get any data in because the number of coordinate pairs in the file is less than the number of cells in the defined region. Markus Metz: Weird. In 6.4, r.in.xyz does import a file where the number of coordinate pairs is far less than the number of cells in the defined region. (I just did a simple test with two input lines and a region with 26.5 million cells, got clean import and correct result) sorry, I'm not really grasping what the problem is. If you suspect there is something weird going on can you make the .xyz.bz2 file available for me to download, and supply the output of r.univar `wc -l` and `r.in.xyz --verbose`? Hamish Ok, the text file (6.1MB compressed) that I use as input to r.in.xyz is available at: http://www.nedbib.za.net/dems/ Output of r.in.xyz --verbose: - Reading data ... Writing to map ... r.in.xyz complete. 0 points found in region. - Output of r.univar: - total null and non-null cells: 1146474 total null cells: 1146474 Of the non-null cells: -- n: 0 minimum: -nan maximum: -nan range: -nan mean: -nan mean of absolute values: -nan standard deviation: -nan variance: -nan variation coefficient: -nan % sum: 0 - output of wc -l: 1146370 Output of r.info: - ++ | Layer:dem_2728ab_25m Date: Mon Jun 7 12:09:34 2010 | Mapset: PERMANENT Login of Creator: hanlie | Location: sa_tm_29deg_E | DataBase: /media/0847147784/grassdata | Title:Raw x,y,z data binned into a raster grid by cell mean ( dem_2728 | Timestamp: none | | | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 1011 | Columns: 1134 | Total Cells: 1146474 |Projection: Transverse Mercator |N: -49312.5S: -74587.5 Res:25 |E: -2987512.5W: -3015862.5 Res:25 | Range of data:min = -nan max = -nan | | Data Source: |/media/0847147784/data/CD-NGI/DEMS/25m dems/2728/ab/2728AB.ort.grass | | | Data Description: |generated by r.in.xyz | | Comments: |r.in.xyz input=/media/0847147784/data/CD-NGI/DEMS/25m dems/2728/ab/\ |2728AB.ort.grass output=dem_2728ab_25m method=mean type=FCELL\ | x=1 y=2 z=3 zscale=1.0 percent=100 - Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: Importing Dems with r.in.xyz
2010/6/7, Markus Metz markus.metz.gisw...@googlemail.com: g.region -p north: -49312.5 south: -74587.5 west: -3015862.5 east: -2987512.5 nsres: 25 ewres: 25 rows: 1011 cols: 1134 cells: 1146474 From the region settings you posted previously, I assume that the 1. column is north and the 2. is east, so: r.in.xyz input=2728AB.ort.grass output=2728AB.ort fs=, x=2 y=1 r.in.xyz complete. 1146370 points found in region. r.univar map=2728AB.ort total null and non-null cells: 1146474 total null cells: 104 Of the non-null cells: -- n: 1146370 minimum: 1473.52 maximum: 1618.85 range: 145.33 mean: 1519.76 mean of absolute values: 1519.76 standard deviation: 26.6106 variance: 708.125 variation coefficient: 1.75098 % sum: 1742202717.5849609375 Everything looks ok. Markus M Ahh, you have found the error. In fact, the columns should not be swopped, the region specifications should. I've retried the import and it works now despite the missing data and the DEM sits in the right place compared to the rest of my data. Thanks for everyone's help and my apologies for the confusion. Our national projection swops X and Y and reverses the direction of the positive coordinates, so it can get very confusing. Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: Importing Dems with r.in.xyz
2010/6/4, Micha Silver mi...@arava.co.il: On 06/04/2010 04:20 PM, Hanlie Pretorius wrote: Hi, I've been following the procedure at http://grass.osgeo.org/wiki/Import_XYZ to import DEMs into GRASS. It has worked for two DEMS, but I have a problem with the third one at the step where one verifies that the number of rows in the ASCII file corresponds to the number of cells in the enlarged region. At this point, g.region reports 1146474 cells in the region, while I have 1146370 lines of coordinates in my file. Assuming you already did the r.in.xyz -s ... step to get and set your region to match the input data. So it looks like there are about 100 coordinates missing from the ASCII file. Maybe holes in the data? One way to work around this might be to import the point data as a 3D vector, then run the usual v.surf.rst interpolation. v.in.xyz -z in=ascii.txt out=elev_pts z=3 fs=, then g.region vect=elev_pts res=choose appropriate resolution v.surf.rst elev_pts elev=dem layer=0 Thanks for the suggestion, Micha. I've tried the interpolation route before and found that the difference between the resulting interpolated surface and the original DEM was up to 7m. I want to use the DEMs for a flood application, so they need to be as accurate as possible. I was thinking perhaps importing the points as vectors, converting them to raster and then doing a nearest neighbour or IDW interpolation to fill the gaps. At least then I'll be able to see where the gaps are and limit the interpolated pixels using a mask? The output of g.region is: - projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -49312.49 south: -74587.5 west: -3015862.5 east: -2987512.5 nsres: 25.0989 ewres: 25 rows: 1011 cols: 1134 cells: 1146474 - The first three and last three coordinates in my text file are: - -74575.00,-3015850.00,1548.83 -74575.00,-3015825.00,1548.33 -74575.00,-3015800.00,1547.50 . . . -49324.99,-2987575.01,1510.98 -49324.99,-2987550.01,1511.24 -49324.99,-2987525.01,1511.47 - Can someone help me to figure out what's going on? 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://www.surfaces.co.il ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Projecting a GCS onto a computer screen
Hi, I was just wondering how the GRASS code displays coordinates in a geographic coordinate system on a computer screen. Do they get projected in some way? Or is there another trick to it? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Problem importing netCDF
Hi, I'm using Win XP and GRASS 6.4RC6 to import a TRMM netCDF file using the command: - r.in.gdal -o input =c:/data/3B42.000202.0.6.nc output=test band=precipitation - and am getting the error: - G_set_window(): Illegal latitude for North - gdalinfo for the file outputs: - Driver: netCDF/Network Common Data Format Files: C:\data\3B42.000202.0.6.nc C:\data\3B42.000202.0.6.nc.aux.xml Size is 512, 512 Coordinate System is `' Metadata: NC_GLOBAL#Conventions=CF-1.0 Subdatasets: SUBDATASET_1_NAME=NETCDF:C:\data\3B42.000202.0.6.nc:precipitation SUBDATASET_1_DESC=[400x1440] precipitation (32-bit floating-point) SUBDATASET_2_NAME=NETCDF:C:\data\3B42.000202.0.6.nc:relativeError SUBDATASET_2_DESC=[400x1440] relativeError (32-bit floating-point) Corner Coordinates: Upper Left (0.0,0.0) Lower Left (0.0, 512.0) Upper Right ( 512.0,0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0) - In Panoply (http://www.giss.nasa.gov/tools/panoply/), the information for the file displays as: - netcdf file:/C:/data/3B42.000202.0.6.nc { dimensions: lat = 400; lon = 1440; variables: float precipitation(lat=400, lon=1440); :long_name = precipitation; :units = mm hr-1; :_FillValue = -.9f; // float float relativeError(lat=400, lon=1440); :long_name = relativeError; :units = mm hr-1; :_FillValue = -.9f; // float float lat(lat=400); :standard_name = latitude; :units = degrees_north; :long_name = latitude; :_CoordinateAxisType = Lat; float lon(lon=1440); :standard_name = longitude; :units = degrees_east; :long_name = longitude; :_CoordinateAxisType = Lon; :Conventions = CF-1.0; } - In the Panoply array display, the coordinates range from 49.875° to -49.875° (lattitude) and -179.875° to +179.875° longitude. These are coordinates for the centres of 0.25° grid cells. I have set my region according to the documentation for the dataset to: - projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 50N south: 50S west: 180W east: 180E nsres: 0:15 ewres: 0:15 rows: 400 cols: 1440 cells: 576000 - Can someone perhaps spot my error? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Importing Dems with r.in.xyz
Hi, I've been following the procedure at http://grass.osgeo.org/wiki/Import_XYZ to import DEMs into GRASS. It has worked for two DEMS, but I have a problem with the third one at the step where one verifies that the number of rows in the ASCII file corresponds to the number of cells in the enlarged region. At this point, g.region reports 1146474 cells in the region, while I have 1146370 lines of coordinates in my file. The output of g.region is: - projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: -49312.49 south: -74587.5 west: -3015862.5 east: -2987512.5 nsres: 25.0989 ewres: 25 rows: 1011 cols: 1134 cells: 1146474 - The first three and last three coordinates in my text file are: - -74575.00,-3015850.00,1548.83 -74575.00,-3015825.00,1548.33 -74575.00,-3015800.00,1547.50 . . . -49324.99,-2987575.01,1510.98 -49324.99,-2987550.01,1511.24 -49324.99,-2987525.01,1511.47 - Can someone help me to figure out what's going on? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: Using an external hard drive as data source
Thanks for all the replies. I reformatted the external drive as FAT and let Ububtu mount it itself, without editing fstab. Then I started GRASS without the sudo option and by selecting the data folder using the /media/volume_name/grassdata path. 2010/5/24, Glynn Clements gl...@gclements.plus.com: Hanlie Pretorius wrote: I'm using GRASS on two machines in different locations (one Win XP, the other Ubuntu 10.04), so I would like to keep my grassdata folder on an external drive. However, when I try to run grass in Ubuntu using the external hard drive (NTFS file system), I get the following error as soon as I've selected the new data location and the GUI starts to load: - Execution failed: 'g.region -u -g -p -c' Details: Error: MAPSET PERMANENT - permission denied - To run GRASS, I use the command 'sudo grass -wx'. GRASS works fine when I run it from the internal hard drive. Is there a way to get GRASS to work from the external drive? You must be the owner of the directory for the current mapset. If you're mounting a FAT filesystem, you can use the uid= mount option to set the ownership of the files. However, I don't know what the situation is with NTFS, as it has its own permission model. -- 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] r.sim.water depth and error ranges are the same
Hi, I ran r.sim.water on a DEM using the process in the manual for generating the rainfall, manning and infiltration maps. The height in the DEM vary from 1429.79 m to 1740.20 m. I opted to output the rainfall depth and the error maps. The ranges of both these output maps (as reported by r.info) are the same: min = 0.000308 max = 0.155246. Is this because of the fake rainfall, manning and infiltration maps, or am I doing something else wrong? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Using an external hard drive as data source
Hi, I'm using GRASS 6.4 on two machines in different locations (one Win XP, the other Ubuntu 10.04), so I would like to keep my grassdata folder on an external drive. However, when I try to run grass in Ubuntu using the external hard drive (NTFS file system), I get the following error as soon as I've selected the new data location and the GUI starts to load: - Execution failed: 'g.region -u -g -p -c' Details: Error: MAPSET PERMANENT - permission denied - To run GRASS, I use the command 'sudo grass -wx'. GRASS works fine when I run it from the internal hard drive. Is there a way to get GRASS to work from the external drive? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Datum unknown in v.in.ogr
Hi, I'm using QGIS to export shapefiles from WGS84 GCS to a local projected coordinate system so that I can import the projected shapefile into GRASS. The PROJ4 definition that I use to export the shapefile (and that I used to create the GRASS location) is: +proj=tmerc +lat_0=0 +lon_0=29 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs When I use v.in.ogr to get the shapfile into GRASS, I get the message: 'Datum unknown not recognised by GRASS and no parameters found Projection of input dataset and current location appear to match' Should I worry about the first part of the message? The imported vector data does look a bit skew. It consists of rectangles that should match raster DEM rectangles, but the outlines now seem to be leaning a few degrees towards the left. Is this just a result of the errors inherent in a projected coordinate system? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: dems from coordinate lists
Hi Jamie, Thanks, this worked perfectly. Regards Hanlie 2010/5/13, Jamie Adams jaad...@gmail.com: Did you adjust the region also? To do this properly, you need to first scan the file using r.in.xyz to get the extents. I'd also use the shell style output. r.in.xyz -sg input=2628cc.ORT.xyz output=dem_2628cc_25m_xyz Should return something like: n=1000 s=500 w=500 e=1000 t=### b=### (you can ignore the t b output) This is the extent of the points, which are going to be the center of the output pixels. Since the region is defined as pixel edge, you need to expand the region 1/2 a pixel in each direction to ensure the points are in the center of the pixel. At 25m spacing, your region should be set as: g.region n=1012.5 s=487.5 w=487.5 e=1012.5 res=25 - Jamie On Thu, May 13, 2010 at 9:02 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi Jamie, I tried r.in.xyz: (r.in.xyz input=2628cc.ORT.xyz output=dem_2628cc_25m_xyz method=mean type=FCELL x=1 y=2 z=3 zscale=1.0 percent=100) but I still get the rows of no data in the result. Hanlie You want the points to represent cell centers so you need to expand your region 1/2 pixel in all four directions. Also look at using r.in.xyz, it works directly on this type of data. Jamie On Thu, May 13, 2010 at 7:26 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I've obtained DEMs in text files with columns X, Y and Z at 25m spacing. The first three entries in the text file are: - X,Y,Z 99550,2.9883e+06,1473.47 99550,2.98828e+06,1473.57 99550,2.98825e+06,1473.63 - To me it seems the easiest way to import these is the following: 1. v.in.ascii 2. Set the region to the extents of new vector file and the resolution to 25m. 3. v.to.rast. This 'works' but I get horizontal strips of no data in my raster DEM at 25m spacings. Also, when I look at the raster and the vector layer together, the vector points are not always on the edges of the raster cells. Can someone perhaps help me to fix this? g.region: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: 2988300 south: 2959850 west: 74300 east: 99550 nsres: 25 ewres: 25 rows: 1138 cols: 1010 cells: 1149380 v.info: ++ | Layer: dem_2628cc_...@c83 | Mapset: C83 | Location:sa_tm_19deg_E | Database:C:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: Administrator | Organization: | Source date: Thu May 13 16:06:07 2010 || | Type of Map: vector (level: 2) | | Number of points: 1151529 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 0 | | Projection: Transverse Mercator | N: 2988300S: 2959850 | E: 99550W: 74300 | B: 1429.79T:1740.2 | | Digitization threshold: 0 | Comments: | ++ r.info: ++ | Layer:dem_2628cc_...@c83 Date: Thu May 13 16:15:50 2010 | Mapset: C83Login of Creator: Administrator | Location: sa_tm_19deg_E | DataBase: C:\Hanlie\grassdata | Title:Categories ( dem_2628cc_25m ) | Timestamp: none || | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1138 | Columns: 1010 | Total Cells: 1149380 |Projection: Transverse Mercator |N:2988300S:2959850 Res:25 |E: 99550W: 74300 Res:25 | Range of data:min = 1429.85 max = 1739.41 | | Data Source: |Vector Map: dem_2628cc_...@c83 in mapset C83 |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=dem_2628cc_...@c83 output=dem_2628cc_25m use=z\ | type=point,line,area layer=1 value=1 rows=4096
Re: [GRASS-user] dems from coordinate lists
Hi Micha, I tried your suggestion after setting the region to 20m instead of the raster DEM's 25m.: v.surf.rst input=dem_2628cc_...@c83 layer=0 elev=dem_2628cc_rst_elev tension=40. segmax=40 npmin=120 dmin=9.998022 dmax=49.990111 zm ult=1.0 This worked, but the differences between the raster DEM that I created with r.in.xyz and the rst interpolated results are quite big - ranging from -6.882202m to +7.864258m. It also ran fairly slowly. Without adjusting the npmin paramter from the default (300) to 120 it literally ran for hours (Win XP, 3GHz CPU, 1GB RAM). Adjusting npmin to 120 didn't seem to affect the error range of the outcome much. Is there a reason why I should use r.surf.rst instead of v.surf.rst? Or perhaps I should just import the points with r.in.xyz and leave the DEM in this format for further applications (hydrological modelling)? Regards Hanlie 2010/5/13, Micha Silver mi...@arava.co.il: MS wrote: If I follow correctly, instead of v.to.rast, you need to interpolate a raster DEM from the points. v.surf.rst produces nice results, but there are other interpolation modules as well in the raster category. That's the method I use also. I start with: v.in.ascii -z in=ascii_file z=3 out=vect_pts This creates a 3D vector using the z column as the height values. Now set the desired region: g.region vect=vect_pts res=xxx Choose the raster resolution that suits your needs. If the points in the ascii file are at 25 m spacing, then you probably could interpolate at 10m-20m resolution (or better) with no problems. Then: v.surf.rst in=vect_pts layer=0 elev=dem ... The layer=0 parameter indicates that you're using the 3D vector's z value for elevation. -- Micha ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] dems from coordinate lists
Hi, I've obtained DEMs in text files with columns X, Y and Z at 25m spacing. The first three entries in the text file are: - X,Y,Z 99550,2.9883e+06,1473.47 99550,2.98828e+06,1473.57 99550,2.98825e+06,1473.63 - To me it seems the easiest way to import these is the following: 1. v.in.ascii 2. Set the region to the extents of new vector file and the resolution to 25m. 3. v.to.rast. This 'works' but I get horizontal strips of no data in my raster DEM at 25m spacings. Also, when I look at the raster and the vector layer together, the vector points are not always on the edges of the raster cells. Can someone perhaps help me to fix this? g.region: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: 2988300 south: 2959850 west: 74300 east: 99550 nsres: 25 ewres: 25 rows: 1138 cols: 1010 cells: 1149380 v.info: ++ | Layer: dem_2628cc_...@c83 | Mapset: C83 | Location:sa_tm_19deg_E | Database:C:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: Administrator | Organization: | Source date: Thu May 13 16:06:07 2010 || | Type of Map: vector (level: 2) | | Number of points: 1151529 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 0 | | Projection: Transverse Mercator | N: 2988300S: 2959850 | E: 99550W: 74300 | B: 1429.79T:1740.2 | | Digitization threshold: 0 | Comments: | ++ r.info: ++ | Layer:dem_2628cc_...@c83 Date: Thu May 13 16:15:50 2010 | Mapset: C83Login of Creator: Administrator | Location: sa_tm_19deg_E | DataBase: C:\Hanlie\grassdata | Title:Categories ( dem_2628cc_25m ) | Timestamp: none || || | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1138 | Columns: 1010 | Total Cells: 1149380 |Projection: Transverse Mercator |N:2988300S:2959850 Res:25 |E: 99550W: 74300 Res:25 | Range of data:min = 1429.85 max = 1739.41 | | Data Source: |Vector Map: dem_2628cc_...@c83 in mapset C83 |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=dem_2628cc_...@c83 output=dem_2628cc_25m use=z\ | type=point,line,area layer=1 value=1 rows=4096 | ++ Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: dems from coordinate lists
Hi Jamie, I tried r.in.xyz: (r.in.xyz input=2628cc.ORT.xyz output=dem_2628cc_25m_xyz method=mean type=FCELL x=1 y=2 z=3 zscale=1.0 percent=100) but I still get the rows of no data in the result. Hanlie You want the points to represent cell centers so you need to expand your region 1/2 pixel in all four directions. Also look at using r.in.xyz, it works directly on this type of data. Jamie On Thu, May 13, 2010 at 7:26 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I've obtained DEMs in text files with columns X, Y and Z at 25m spacing. The first three entries in the text file are: - X,Y,Z 99550,2.9883e+06,1473.47 99550,2.98828e+06,1473.57 99550,2.98825e+06,1473.63 - To me it seems the easiest way to import these is the following: 1. v.in.ascii 2. Set the region to the extents of new vector file and the resolution to 25m. 3. v.to.rast. This 'works' but I get horizontal strips of no data in my raster DEM at 25m spacings. Also, when I look at the raster and the vector layer together, the vector points are not always on the edges of the raster cells. Can someone perhaps help me to fix this? g.region: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: 2988300 south: 2959850 west: 74300 east: 99550 nsres: 25 ewres: 25 rows: 1138 cols: 1010 cells: 1149380 v.info: ++ | Layer: dem_2628cc_...@c83 | Mapset: C83 | Location:sa_tm_19deg_E | Database:C:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: Administrator | Organization: | Source date: Thu May 13 16:06:07 2010 || | Type of Map: vector (level: 2) | | Number of points: 1151529 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 0 | | Projection: Transverse Mercator | N: 2988300S: 2959850 | E: 99550W: 74300 | B: 1429.79T:1740.2 | | Digitization threshold: 0 | Comments: | ++ r.info: ++ | Layer:dem_2628cc_...@c83 Date: Thu May 13 16:15:50 2010 | Mapset: C83Login of Creator: Administrator | Location: sa_tm_19deg_E | DataBase: C:\Hanlie\grassdata | Title:Categories ( dem_2628cc_25m ) | Timestamp: none || | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1138 | Columns: 1010 | Total Cells: 1149380 |Projection: Transverse Mercator |N:2988300S:2959850 Res:25 |E: 99550W: 74300 Res:25 | Range of data:min = 1429.85 max = 1739.41 | | Data Source: |Vector Map: dem_2628cc_...@c83 in mapset C83 |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=dem_2628cc_...@c83 output=dem_2628cc_25m use=z\ | type=point,line,area layer=1 value=1 rows=4096 | ++ Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user -- next part -- An HTML attachment was scrubbed... URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20100513/dd33d293/attachment.html -- ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user End of grass-user Digest, Vol 49, Issue 24 ** ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: maps France and Germany
Try http://www.gadm.org. Date: Wed, 12 May 2010 06:23:34 +0200 From: Jos? Miguel Barrios jmbarri...@gmail.com Subject: [GRASS-user] maps France and Germany To: grass-user@lists.osgeo.org Message-ID: aanlktimmcisrcoamohjfuhsbdi2l7wdqejzfsr2ra...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Hi list; Can someone indicate me where can I request simple vector layers of Germany and France at municipality level? Thanks! Miguel -- next part -- An HTML attachment was scrubbed... URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20100512/688ec0da/attachment.html -- ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user End of grass-user Digest, Vol 49, Issue 20 ** ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] d.* command problems
Hi Glynn, This worked, thanks. 2010/5/11, Glynn Clements gl...@gclements.plus.com: Hanlie Pretorius wrote: I just pasted the code you sent at the end of the d.correlate script That won't work. The code is a patch, indicating changes to be made to the file. I've attached a complete modified version of the script. -- Glynn Clements gl...@gclements.plus.com #!/bin/sh # # # MODULE: d.correlate for GRASS 6; based on dcorrelate.sh for GRASS 4,5 # AUTHOR(S):CERL - Michael Shapiro; updated to GRASS 6 by Markus Neteler 5/2005 # PURPOSE: prints a graph of the correlation between data layers (in pairs) # derived from grass5/src.local/d.correlate.sh # COPYRIGHT:(C) 2005 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. # # TODO GRASS 7: rename parameters to map1, map2 or better multiple map=map1[,map2[...]] # #%Module #% description: Prints a graph of the correlation between data layers (in pairs). #% keywords: display, diagram #%End #%option #% key: layer1 #% type: string #% gisprompt: old,cell,raster #% description: raster input map #% required : yes #%end #%option #% key: layer2 #% type: string #% gisprompt: old,cell,raster #% description: raster input map #% required : yes #%end #%option #% key: layer3 #% type: string #% gisprompt: old,cell,raster #% description: raster input map #% required : no #%end #%option #% key: layer4 #% type: string #% gisprompt: old,cell,raster #% description: raster input map #% required : no #%end if [ -z $GISBASE ] ; then echo You must be in GRASS GIS to run this program. 2 exit 1 fi if [ $1 != @ARGS_PARSED@ ] ; then CMDLINE=`basename $0` for arg in $@ ; do CMDLINE=$CMDLINE \$arg\ done export CMDLINE exec g.parser $0 $@ fi PROG=`basename $0` check if we have awk if [ ! -x `which awk` ] ; then g.message -e awk required, please install awk or gawk first exit 1 fi # setting environment, so that awk works properly in all languages unset LC_ALL LC_NUMERIC=C export LC_NUMERIC if [ $# -gt 4 ] ; then g.message -e max 4 layers allowed exit 1 fi TMP1=`g.tempfile pid=$$` ok=yes for f in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 do eval `g.findfile element=cell file=$f` if [ -z $name ] ; then g.message -e $f not found ok=no fi done if [ $ok = no ] ; then exit 1 fi d.erase GRASS_PNG_READ=TRUE export GRASS_PNG_READ if [ $? != 0 ] ; then exit 1 fi # how many cmd line arguments? ARGNUM=`echo $CMDLINE | tr -s ' ' '\n' | wc -l | awk '{print $1 - 1}'` echo CORRELATION | d.text color=white size=4 line=1 colors=red black blue green gray violet if [ $ARGNUM -eq 2 ] ; then line=4 else line=2 fi iloop=0 jloop=0 # get max in case of two maps for x, y axis eval `r.univar -g $GIS_OPT_LAYER1` max_layer1=$max eval `r.univar -g $GIS_OPT_LAYER2` max_layer2=$max for i in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 do iloop=`expr $iloop + 1` for j in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do jloop=`expr $jloop + 1` if [ $i != $j -a $iloop -le $jloop ] ; then colorstmp1=`echo $colors | cut -d' ' -f1` colorstmp2=`echo $colors | cut -d' ' -f2-` colors=`echo $colorstmp2 $colorstmp1` if [ $ARGNUM -eq 2 ] ; then echo $j | d.text color=`echo $colors | cut -d' ' -f1` size=4 at=0,9$line echo $i | d.text color=`echo $colors | cut -d' ' -f1` size=4 at=60,0$line else echo $i $j | d.text color=`echo $colors | cut -d' ' -f1` size=4 line=$line fi line=`expr $line + 1` r.stats -cnA input=$i,$j $TMP1 m=`awk '$1max1{max1=$1} $2max2{max2=$2} min1==0||$1min1{min1=$1} min2==0||$2min2{min2=$2} END {print min1,max1,min2,max2}' $TMP1` m1=`echo $m | cut -d' ' -f1` m2=`echo $m | cut -d' ' -f2` m3=`echo $m | cut -d' ' -f3` m4=`echo $m | cut -d' ' -f4` awk '{print move,($1-min1+1)*100.0/(max1-min1+1),($2-min2+1)*100.0/(max2-min2+1);print draw,($1-min1+1)*100.0/(max1-min1+1),($2-min2+1)*100.0/(max2-min2+1) }' min1=$m1 max1=$m2 min2=$m3 max2=$m4 $TMP1 | d.graph color=`echo $colors | cut -d' ' -f1` if [ $ARGNUM -eq 2 ] ; then d.graph EOF size 2 2 move 0 92 text $max_layer1 move 90 2 text $max_layer2 EOF fi fi done jloop=0 done rm -f $TMP1 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] d.* command problems
Hi, I'm running GRASS 6.4.0RC6 on Windows XP SP1. Working my way through the book Open Source GIS a GRASS GIS Approach 3rd ed. by Markus Neteler and Helena Mitasova, I have noticed problems with some d.* commands. For example, when I type d.histogram into the GUI layer manager command line, I get the response Command 'd.histogram' not yet implemented in the GUI. Try adding it as a command layer instead. When I add the command layer, I can't see what I'm typing and I can't find a way to adjust the display so that I can see what's going on. d.histogram is still ok because I can view the output in the msys folder. However, d.correlate seems to output multiple files that overwrite one another. Is there a way around this? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.contour generates no lines
Hi, I tried to create contours from a DEM with the following information: | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 9856 | Columns: 4844 | Total Cells: 47742464 |Projection: Transverse Mercator |N: -3704977.57894189S: -3803537.57894189 Res:10 |E: -15589.63705055W: -64029.63705055 Res:10 | Range of data:min = -39.681747 max = 1084.584839 The command I used was: r.contour --overwrite input=ct_dem_...@permanent output=ct_contours_1...@permanent step=100 cut=2000 However, the resulting vector file contains no lines. This is the output of v.info: | Type of Map: vector (level: 2) | | Number of points: 0 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 1 | | Projection: Transverse Mercator | N: 0S: 0 | E: 0W: 0 | B: 0T: 0 The attribute table does contain 10 categories starting from 100m and ending at 1000m. Leaving out the cut parameter doesn't change anything. Any ideas about why the file contains no features? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Problem after import using r.in.ascii
Greetings Herbivores, This is the first time that I've created a location and mapset fo my own and I'm trying to import TRMM precipitation data into GRASS using r.in.ascii. I have an existing CSV file that has the headers lat, long and precip. The precipitation is a floating point number to two decimal places. The coordinates in the CSV file define the centre of a quarter degree square. The first set of coordinates in the file is -19.625,15.875. Becuase the TRMM data defines the centre of the square, I subtracted 0.125 from the first set of coordinates to get the north-west corner of the raster. And I added 0.125 to the last coordinate in the CSV file, which is -34.875,34.125, to get the south-east corner of the raster. In the CSV file, the latitude stays constant for 74 rows while the longitude increments in quarter degree steps. The whole file contains 4588 rows, so a grid of the data should have 74 rows and 62 columns, right? So I reformatted the data to 74 rows and 62 columns and added the header: north: -19.5 south: -35.0 east:15.75 west:34.25 rows:74 cols:62 null:-319.99 The output of g.region -p is: projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 18S south: 36S west: 15E east: 35E nsres: 0:15 ewres: 0:15 rows: 72 cols: 80 cells: 5760 I then imported the reformatted file using the command: r.in.ascii -f input=3B42RT.2010032409.6.bin.sa.grass\ output=2010032409.6.precip \ title=Precipitation 2010032409.6\ mult=1.0 or read from header nv=* or read from header When imported, I get two thin strips of data on the sides of the image. The raster shouldn't contain any nulls (the text file doesn't contain the value -319.99), so something is wrong. r.info gives a very strange resolutions: | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 74 | Columns: 62 | Total Cells: 4588 |Projection: Latitude-Longitude |N: 19:30SS:35S Res: 0:12:34.054054 |E: 15:45EW: 34:15E Res: 5:30:29.032258 | Range of data:min = 0.00 max = 7.47 Can someone perhaps help me to see what's happening here? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.sim.water crashes (repost)
Hi, This is a repost since I got no reaction the previous time. I'm not sure if it is because people can't help. Is it ok to contact the authors at the bottom of the r.sim.water page if no one on this list can help me? I'm working on Windows XP SP3 and Grass 6.4.0RC6. When I try to run r.sim.water according to the example on page 163 in the book Open Source GIS a GRASS GIS approach (3rd Ed.), I get one of those Windows errors that want to 'tell Microsoft about this problem'. The details of the error report says: AppName: r.sim.water.exe AppVer: 0.0.0.0 ModName: libgrass_gis.6.4.0rc6.dll ModVer: 0.0.0.0 Offset: 0001a3f7 The first module listed in the detailed technical report is r.sim.water.exe. Many modules follow after that. The book uses the North Carolina example datasets and the steps are as follows: g.region=rural_1m res=2 -p r.mapcalc man05=0.05 r.mapcalc infil0=0. r.mapcalc rain50mmhr=50 v.surf.rst -d input=elev_lid792_bepts layer=0 elev=elev_lid792_2m slope=dx_2m aspect=dy_2m ten=15 smooth=1.5 segmax=25 npmin=100 r.sim.water -t elevin=elev_lid792_2m dxin=dx_2m dyin=dy_2m rain=rain50mmhr infil=infil0 manin=man05 depth=wdp_2m disch=disch_2m nwalk=40 niter=500 outit=20 hmax=0.2 halpha=8.0 hbeta=1.0 I have tried these steps on two different computers and get the same error. I have also tried the example on the manual page of r.sim.water with the same result. Can someone perhaps suggest a fix? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: r.sim.water crashes (repost)
I suppose this must be related to ticket 617 (http://trac.osgeo.org/grass/ticket/617). Guess I'll have to switch to Linux if I want to work with it now. 2010/4/9, Hanlie Pretorius hanlie.pretor...@gmail.com: Hi, This is a repost since I got no reaction the previous time. I'm not sure if it is because people can't help. Is it ok to contact the authors at the bottom of the r.sim.water page if no one on this list can help me? I'm working on Windows XP SP3 and Grass 6.4.0RC6. When I try to run r.sim.water according to the example on page 163 in the book Open Source GIS a GRASS GIS approach (3rd Ed.), I get one of those Windows errors that want to 'tell Microsoft about this problem'. The details of the error report says: AppName: r.sim.water.exe AppVer: 0.0.0.0 ModName: libgrass_gis.6.4.0rc6.dll ModVer: 0.0.0.0 Offset: 0001a3f7 The first module listed in the detailed technical report is r.sim.water.exe. Many modules follow after that. The book uses the North Carolina example datasets and the steps are as follows: g.region=rural_1m res=2 -p r.mapcalc man05=0.05 r.mapcalc infil0=0. r.mapcalc rain50mmhr=50 v.surf.rst -d input=elev_lid792_bepts layer=0 elev=elev_lid792_2m slope=dx_2m aspect=dy_2m ten=15 smooth=1.5 segmax=25 npmin=100 r.sim.water -t elevin=elev_lid792_2m dxin=dx_2m dyin=dy_2m rain=rain50mmhr infil=infil0 manin=man05 depth=wdp_2m disch=disch_2m nwalk=40 niter=500 outit=20 hmax=0.2 halpha=8.0 hbeta=1.0 I have tried these steps on two different computers and get the same error. I have also tried the example on the manual page of r.sim.water with the same result. Can someone perhaps suggest a fix? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.sim.water crashes
Hi, I'm working on Windows XP SP3 and Grass 6.4.0RC6. When I try to run r.sim.water according to the example on page 163 in the book Open Source GIS a GRASS GIS approach (3rd Ed.), I get one of those Windows errors that want to 'tell Microsoft about this problem'. The details of the error report says: AppName: r.sim.water.exe AppVer: 0.0.0.0 ModName: libgrass_gis.6.4.0rc6.dll ModVer: 0.0.0.0 Offset: 0001a3f7 The first module listed in the detailed technical report is r.sim.water.exe. Many modules follow after that. The book uses the North Carolina example datasets and the steps are as follows: g.region=rural_1m res=2 -p r.mapcalc man05=0.05 r.mapcalc infil0=0. r.mapcalc rain50mmhr=50 v.surf.rst -d input=elev_lid792_bepts layer=0 elev=elev_lid792_2m slope=dx_2m aspect=dy_2m ten=15 smooth=1.5 segmax=25 npmin=100 r.sim.water -t elevin=elev_lid792_2m dxin=dx_2m dyin=dy_2m rain=rain50mmhr infil=infil0 manin=man05 depth=wdp_2m disch=disch_2m nwalk=40 niter=500 outit=20 hmax=0.2 halpha=8.0 hbeta=1.0 I have tried these steps on two different computers and get the same error. I have also tried the example on the manual page of r.sim.water with the same result. Can someone perhaps suggest a fix? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Language in GRASS help files
Hi, I've just installed GRASS 6.4.0svn for Windows (Running on Windows XP SP3) and I find a peculiar thing happening in the help files. I was looking at the help page for the TOPMODEL module and at the top of the page the text displays in a Slavic language. Further down on the page some English text displays as well. The location on the computer is set to English (South Africa). Does someone perhaps know how I can get all the help text in English? Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user