Micha, I followed your instructions with the following results:
----
Import data from 2 shapefiles, Contour (layer 1)and Track (layer 2): (no
problems reported)
v.in.ogr dsn=/home/george/GRASSDATA/ output=testvectmap min_area=0.0001
snap=-1
----
Display vector map: (no problems, map displayes correctly)
d.vect map=testvectmap@testvector color=0:0:0 lcolor=0:0:0
fcolor=170:170:170 display=shape type=point,line,boundary,area
icon=basic/circle size=5 layer=-1 lsize=8 xref=left yref=center llayer=1
----
Convert vector to raster and thin: (no problems reported)
v.to.rast input=testvectmap@testvector output=testrastmap use=attr
type=point,line,area layer=1 column=ELEVATION value=1 rows=100000
r.thin input=testrastmap@testvector output=thinrastmap iterations=200
----
Display raster map (PROBLEM:only displays large, brightly coloured
rectangles, no contours
d.rast map=testrastmap@testvector -o
----
At this point I stopped as I got the impression something was already wrong
but could not work out what.
Also, should I not rasterize the track data as well?
Any help would be appreciated. If you require the source shapefiles I have
no problems sending them to you.
Micha Silver wrote:
Yes.
I think you should be looking at converting the shapefile contour lines
to a raster....
v.in.org to import the contour lines and tracks into grass
v.to.rast and r.thin to convert the contour lines to a raster map.
r.surf.contour to create an elevation surface from the rasterized
contours
r.profile to get the elevations along the tracks.
Post back with specific questions if you hit a snag.
Micha, I followed your instructions with the following results:
----
Import data from 2 shapefiles, Contour (layer 1)and Track (layer 2): (no
problems reported)
v.in.ogr dsn=/home/george/GRASSDATA/ output=testvectmap min_area=0.0001
snap=-1
----
How many shapefiles did you have in the GRASSDATA directory?
What does testvectmap contain? You should make two separate GRASS vectors, one for the tracks and one for the contours. You can use the ‘layer’ parameter each time to specify each shapefile, like:
v.in.ogr dsn=/home/george/GRASSDATA/ layer=contours.shp output=contours
At this point I’d check to make sure the elevation data came thru:
v.info -c # shows the column headers, to make sure you have ELEVATION
v.db.select # shows the actual ELEVATION values
Convert vector to raster and thin: (no problems reported)
v.to.rast input=testvectmap@testvector output=testrastmap use=attr
type=point,line,area layer=1 column=ELEVATION value=1 rows=100000
r.thin input=testrastmap@testvector output=thinrastmap iterations=200
----
Display raster map (PROBLEM:only displays large, brightly coloured
rectangles, no contours
d.rast map=testrastmap@testvector -o
----
How did you set up the Location/Mapset ? You might have a problem with the resolution.
What do you get from:
r.info testrastmap
g.region -p
At this point I stopped as I got the impression something was already wrong
but could not work out what.
Also, should I not rasterize the track data as well?
Any help would be appreciated. If you require the source shapefiles I have
no problems sending them to you.
Let’s keep stabbing at it from your end. That’s always the best way to understand the methodology.
Micha Silver wrote:
How many shapefiles did you have in the GRASSDATA directory? …
I had the two shapefiles, Contour and Track, but I got your message loud and clear about creating separate vectors, which I did and displayed as follows: v.in.ogr -o dsn=/home/george/GRASSDATA/Contour.shp output=contours min_area=0.0001 type=line snap=-1 --overwrite v.in.ogr -o dsn=/home/george/GRASSDATA/Track_cl.shp output=tracks min_area=0.0001 type=line snap=-1 --overwrite d.vect map=contours@test02 color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=point,line,boundary,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1 d.vect map=tracks@test02 color=255:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=point,line,boundary,area icon=basic/circle size=5 width=3 layer=1 lsize=8 xref=left yref=center llayer=1 No problems so far, all displayed correctly.
Micha Silver wrote:
At this point I’d check to make sure the elevation data came thru:
v.info -c contours@test02 Displaying column types/names for database connection of layer 1: INTEGER|cat DOUBLE PRECISION|ID CHARACTER|CRT_DATE CHARACTER|MOD_DATE DOUBLE PRECISION|ELEVATION CHARACTER|DESIGNATED CHARACTER|NAT_FORM CHARACTER|DEFINITION v.info -c tracks@test02 Displaying column types/names for database connection of layer 1: INTEGER|cat DOUBLE PRECISION|ID CHARACTER|NAME CHARACTER|CRT_DATE CHARACTER|MOD_DATE CHARACTER|TRACK_USE CHARACTER|TRACK_TYPE CHARACTER|TRK_STATUS DOUBLE PRECISION|area DOUBLE PRECISION|len v.db.select map=contours@test02 layer=1 column=ELEVATION fs=| ELEVATION 1180 1220 1200 etc… etc… So far all is as it should be (I think…) Convert vector to raster: (no problems reported) v.to.rast input=contours@test02 output=contours_rast use=attr type=line layer=1 column=ELEVATION value=1 rows=100000 --overwrite (… Converted points/lines: 630 of 630) ---- Display raster map (PROBLEM:only displays large, brightly coloured rectangles, no contours d.rast map=contour_rast@test02 -o ----
Hamish, please forgive my ignorance but I am not sure what you mean. As far
as I know g.region is the means to change the resolution, whether through
the command line or the GUI. Can you please explain how to access and change
the computational region settings.
hamish_b wrote:
[please don't post in html, it makes it hard to read]
georgew wrote:
PS I tried to play around with resolution in the Region dialog but the
display never changed from the same chunky coloured squares.
you need to change the computational region, not just the display region
settings. the computational region is what the processing modules will
use.
Hamish, please forgive my ignorance but I am not sure what you mean.
As far as I know g.region is the means to change the resolution,
correct,
whether through the command line or the GUI. Can you please explain how
to access and change the computational region settings.
I was talking about when using it from the GUI.
from the command line the display (d.mon) resolution is the computation
region.
in the GUI there are two: the active display region(+resolution) and the
computational one. What you see in the display window(s) is independent
of what's in the $MAPSET/WIND file. That WIND(ow) file holds the region
which will be used by any modules running. The reason for all this
confusion is that you can have multiple GUI display windows open at the
same time, each independently zoomed to a different area and set at a
different resolution (for a visual display there's no point at being at a
higher resolution than you can see). Multiple displays would need multiple
WIND files but the processing module (eg v.to.rast) wouldn't know which
was the appropriate one to use. Thus computational region refers to what
is stored in the mapset's WIND file.
(sorry gui programmers if I have made small mistakes in that)
In gis.m, from the Map Display window if you click on the magnifying glass
icon (with the map) you get a pulldown menu. the last item is "set computational region extents to match display".
also in the GIS manager menu there is Config->Region. In that is "Display
region settings"* and "Change region settings" which brings up the
g.region GUI as if you were running it at the command prompt.
[*] perhaps badly worded, nothing to do with the Map Display- it prints
the current computational region to the output window.
A symptom of the computational region being way out of whack compared to
your display region is that you see the result as a series of huge blocks
(computational resolution is too coarse), or in the other direction the
resolution is so fine that takes ages to load. YMMV
Hamish, please forgive my ignorance but I am not sure what you mean.
As far as I know g.region is the means to change the resolution,
correct,
whether through the command line or the GUI. Can you please explain how
to access and change the computational region settings.
I was talking about when using it from the GUI.
from the command line the display (d.mon) resolution is the computation
region.
in the GUI there are two: the active display region(+resolution) and the
computational one. What you see in the display window(s) is independent
of what's in the $MAPSET/WIND file. That WIND(ow) file holds the region
which will be used by any modules running. The reason for all this
confusion is that you can have multiple GUI display windows open at the
same time, each independently zoomed to a different area and set at a
different resolution (for a visual display there's no point at being at a
higher resolution than you can see). Multiple displays would need multiple
WIND files but the processing module (eg v.to.rast) wouldn't know which
was the appropriate one to use. Thus computational region refers to what
is stored in the mapset's WIND file.
(sorry gui programmers if I have made small mistakes in that)
In gis.m, from the Map Display window if you click on the magnifying glass
icon (with the map) you get a pulldown menu. the last item is "set computational region extents to match display".
And, the reverse: Zoom display to computational region. Maybe your v.to.rast worked as it should, but you just don't see the correct region in the Map Display ?
BTW:
I set the Location in GRASS's opening welcome screen, using the Contour.shp file to provide the boundaries.
and
g.region -p projection: 0 (x,y)
You are aware that this means that the Contour.shp did not contain any projection info and that GRASS thus consideres it as unprojected ?
A symptom of the computational region being way out of whack compared to
your display region is that you see the result as a series of huge blocks
(computational resolution is too coarse), or in the other direction the
resolution is so fine that takes ages to load. YMMV
Note that a mismatch between computational and display regions isn't
necessarily "wrong". You may want to perform computations at a
resolution finer than that of the display for accuracy, or coarser
than it for speed.
Micha Silver wrote:
How many shapefiles did you have in the GRASSDATA directory? ....
provide the boundaries. Here is the result: r.info contours_rast +----------------------------------------------------------------------------+ | Layer: contours_rast Date: Thu Dec 4 08:50:19 2008 | | Mapset: test02 Login of Creator: george | | Location: test02 | | DataBase: /home/george/GRASSDATA | | Title: Labels ( contours_rast ) | | Timestamp: none | |----------------------------------------------------------------------------| | | | Type of Map: raster Number of Categories: 0 | | Data Type: DCELL | | Rows: 20 | | Columns: 20 | | Total Cells: 400 | | Projection: x,y | | N: 6020724.5 S: 6013982.586 Res: 337.0957 | | E: 2481099.136 W: 2467034.878 Res: 703.2129 | | Range of data: min = 520.000000 max = 1840.000000 | | | | Data Source: | | Vector Map: contours@test02 in mapset test02 | | Original scale from vector map: 1:1 | | | | Data Description: | | generated by v.to.rast | | | | Comments: | | v.to.rast input="contours@test02" output="contours_rast" use="attr" \ | | type="line" layer=1 column="ELEVATION" value=1 rows=100000 | | | +----------------------------------------------------------------------------+ GRASS 6.3.0 (test02):~ > g.region -p projection: 0 (x,y) zone: 0 north: 6020724.5 south: 6013982.586 west: 2467034.878 east: 2481099.136 nsres: 337.0957 ewres: 703.2129 rows: 20 cols: 20 cells: 400 I hope all this helps. Your help so far has been invaluable, thanks a lot. PS I tried to play around with resolution in the Region dialog but the display never changed from the same chunky coloured squares.
Here's What I can make out from the above:
The X-Y values are about 6,000,000 N and about 2,450,000 E, and the size of your region is about 14,000 X 6500. Let's *assume* these numbers are in meters (UTM projection??). But your resolution is set to 337 e-w and 703 n-s, thus giving you a raster of 20X20 cells. That's certainly too coarse.
I'd suggest:
Blow away the rasters you made so far.
Set resolution to 10m. X 10m. with
g.region -s -p res=10
You should then see about 1400 columns by 650 rows =~ 910000 cells.
Now recreate those rasters. You should get better looking results.
Micha Silver wrote:
Set resolution to 10m. X 10m. with g.region -s -p res=10
Hamish_b wrote:
Multiple displays would need multiple WIND files but the processing module (eg v.to.rast) wouldn’t know which was the appropriate one to use. Thus computational region refers to what is stored in the mapset’s WIND file.
Moritz Lennert-2 wrote:
You are aware that this means that the Contour.shp did not contain any projection info and that GRASS thus consideres it as unprojected ?
Setting the resolution BEFORE converting to raster was the answer, thank you. I now see the contours in the raster map as well as in the original vector file. As for the null projection, thanks for pointing it out (I was not aware!) but since the purpose of the exercise was to create an elevation profile from tabular data perhaps it is not critical. So I went ahead with the original purpose, i.e. to create an elevation profile for the track.
Micha Silver wrote:
r.surf.contour to create an elevation surface from the rasterized contours
r.profile to get the elevations along the tracks.
So I entered:
r.surf.contour input=contours_rast@PERMANENT output=elevmap_raster --overwrite
and few seconds later it completed its job without errors. Then went ahead with:
r.profile -i input=elevmap_rast@PERMANENT output=- null=*
I followed the track with the left mouse button and a two column list of values was printed on stdout, the distance and the elevation. Success!
My only remaining question is: Is there a way to use the track raster file to input the values instead of using the mouse, so that the process can be more or less automated? Thank you all for your great help. I’ll write a summary and post it in the Wiki as an FAQ or short tutorial as suggested by Markus Neteler.
Much progress!
My only remaining question is: Is there a way to use the track raster file to input the values instead of using the mouse, so that the process can be more or less automated? Thank you all for your great help. I'll write a summary and post it in the Wiki as an FAQ or short tutorial as suggested by Markus Neteler.
Reading the r.profile man page, you will see that it's possible to pipe the list of coordinates to r.profile:
The *profile* parameter can be set to comma separated geographic coordinates for profile line endpoints. The interactive flag (*-i*) overrides this option. Alternatively the coordinate pairs can be piped from stdin, one comma separated pair per line.
You can get the coordinate pairs from your track vector with v.out.ascii, and use the fs=, option to make the output file "comma-separated".
Much progress!
My only remaining question is: Is there a way to use the track raster file
to input the values instead of using the mouse, so that the process can be
more or less automated? Thank you all for your great help. I'll write a
summary and post it in the Wiki as an FAQ or short tutorial as suggested by
Markus Neteler.
I've used the command "v.drape" to get profiles on raster maps along
lines from vector maps. The command to be executed is something like.