[GRASS-user] import satellite data in Grass

Martina_Schaefer · June 4, 2008, 6:12pm

Hi,

I'm new with Grass and I don't have somebody around me knowing well Grass and besides all the documentation and tutorials that I've found, I have some beginner questions...

I have several satellite data sets, that I would like to import in Grass.

These sets are:
- devpan_2000_mos.tif, .rrd, .aux, .tif.aux
- S4_07935_7537_20070828_p10_1_lcc00.tif, .xml with files selected_zone-zone_selectionnee.dbf, prj, shp and shx as well as
spot_footprint-empreinte_spot.dbf, prj, shp and shx
- oa16787_087.img

What I am doing is:
- enter Grass with the spearfish60 tutorials
- then trying to create a new location
in the first case this is:
r.in.gdal -e devpan_2000_mos.tif location=DevonNew output=testfile
and I get this error
ERROR 4: `devpan_2000_mos.tif.aux' not recognised as a supported file format.
A datum name wgs84 (WGS_1984) was specified without transformation parameters.
Note that the GRASS default for wgs84 is towgs84=0.000,0.000,0.000.
100%
CREATING SUPPORT FILES FOR testfile
SETTING GREY COLOR TABLE FOR testfile (8bit, full range)
r.in.gdal complete.

in the second case this is:
  r.in.gdal -e ../Geobase/s4_07935_7537_20070828_p10_lcc00/S4_07935_7537_20070828_p10_1_lcc00.tif location=DevonGeo output=testfile2
and I get once more somme warnings
WARNING: Datum 'Not_specified_based_on_GRS_1980_ellipsoid' not recognised
          by GRASS and no parameters found. Datum transformation will not be
          possible using this projection information.
  100%
CREATING SUPPORT FILES FOR testfile2
SETTING GREY COLOR TABLE FOR testfile2 (8bit, full range)
r.in.gdal complete.

in the third cas this is:
r.in.gdal -e oa16787_087.img location=DevonIMG output=test3
with the following error
ERROR 4: `oa16787_087.img' not recognised as a supported file format.

In the first two cases I can visualize my data then by exiting and starting GRASS again.
But I'm not sure about the warnings and errors that I got. And in the third case I can't access my data.

Can somebody help me with these questions?

And, is there a way to open several files in the same time? In fact, I don't have only these 3 types of data, but also several maps of these type that I would like to put together or at least to visualize together without leaving and starting Grass every time!

I'll have some more questions about the analaysis that I want to do later, but hopefully these problems could be solved first!

Thanks a lot for your help,

Martina

Hamish · June 5, 2008, 1:13am

Martina Schaefer wrote:

I'm new with Grass and I don't have somebody around me knowing well
Grass and besides all the documentation and tutorials that
I've found, I have some beginner questions...

I have several satellite data sets, that I would like to
import in Grass.

These sets are:
- devpan_2000_mos.tif, .rrd, .aux, .tif.aux
- S4_07935_7537_20070828_p10_1_lcc00.tif, .xml with files
selected_zone-zone_selectionnee.dbf, prj, shp and shx as
well as
spot_footprint-empreinte_spot.dbf, prj, shp and shx
- oa16787_087.img

first step: have a look at their metadata with "gdalinfo <filename>"
check that is all in order.

if projections are all the same use r.in.gdal from the new location to load in the other maps. If they use different projections you can use r.proj to move them into the same location. each different map projection needs its own location.

qgis is a nice data viewer if you don't want to create new grass locations just to look at the data. (www.qgis.org)

Hamish

Martina_Schaefer · June 19, 2008, 8:30pm

Hi again,

thanks for all your fast answers to my beginner questions :-).

In the meanwhile I managed to import my data, even the img files, but only by transforming them into geoTIF.
I'm still missing some information about the projections parameters, but for the moment I can use the default values.

What I want to do now, is to overlay the images with a raster and define (in a first time at least by hand) a value for each of the the raster cells.
Is there anything foreseen in Grass that could help me with this? Or do I have just to write a table in any other editor?

What I also want to do, is measure distances and angles between objects on my images. For the moment I don't really see where to search for this kind of tools. I thought that it would be possible to draw lines and measure them and angles between them?

Thanks again for your helpful comments.

Martina

Hamish wrote:

Martina Schaefer wrote:

I'm new with Grass and I don't have somebody around me knowing well Grass and besides all the documentation and tutorials that
I've found, I have some beginner questions...

I have several satellite data sets, that I would like to
import in Grass.

These sets are:
- devpan_2000_mos.tif, .rrd, .aux, .tif.aux
- S4_07935_7537_20070828_p10_1_lcc00.tif, .xml with files selected_zone-zone_selectionnee.dbf, prj, shp and shx as
well as
spot_footprint-empreinte_spot.dbf, prj, shp and shx
- oa16787_087.img

first step: have a look at their metadata with "gdalinfo <filename>"
check that is all in order.

if projections are all the same use r.in.gdal from the new location to load in the other maps. If they use different projections you can use r.proj to move them into the same location. each different map projection needs its own location.

qgis is a nice data viewer if you don't want to create new grass locations just to look at the data. (www.qgis.org)

Hamish

Hamish · June 20, 2008, 1:06pm

Martina:

What I want to do now, is to overlay the images with a raster and
define (in a first time at least by hand) a value for each of the
the raster cells. Is there anything foreseen in Grass that could
help me with this? Or do I have just to write a table in any other
editor?

could you please expand a little bit more? I have trouble understanding exactly the process you wish to do.. r.coin? r.patch? r.mapcalc?
'd.rast -o' + d.what.rast?

What I also want to do, is measure distances and angles between
objects on my images. For the moment I don't really see where
to search for this kind of tools. I thought that it would be possible
to draw lines and measure them and angles between them?

m.cogo might help. or 'v.distance upload=to_along,to_angle', d.measure, v.to.db...

In general beware to double check results/meaning if using a lat/lon location (see d.rhumbline and d.geodesic), or map projection with possibly heavy distortion at far distances.

If you build the Python SWIG interface there is also a m.distance module to look at.

see also this thread+code patch from March:
"How to find out an angle between to points on the map"
http://thread.gmane.org/gmane.comp.gis.grass.user/22817/
http://thread.gmane.org/gmane.comp.gis.grass.devel/25396/

Hamish

Martina_Schaefer · June 20, 2008, 1:33pm

Hi Hamish,

thanks for your fast answer.

Hamish wrote:

Martina:

What I want to do now, is to overlay the images with a raster and
define (in a first time at least by hand) a value for each of the
the raster cells. Is there anything foreseen in Grass that could
help me with this? Or do I have just to write a table in any other
editor?

could you please expand a little bit more? I have trouble understanding exactly the process you wish to do.. r.coin? r.patch? r.mapcalc?
'd.rast -o' + d.what.rast?

What I want to do is the following:
I have satellite photos from a glacier and we want to get an idea of the distribution of the crevasses over this glacier.
In a first step we want to assign every 100x100m grid cell of this glacier a value reflecting if there are very little, some, a lot of crevasses. As we don't know for the moment how to do this automatically (perhaps it would be possible to use the color in a next step), the idea was to define this value by "eyeballing". Any ideas how to do this a better way are welcome, too, but I think it's not easy.
Now the questions is "only" how to handle the grid and my values in Grass ... or perhaps it would be better to print the image + the grid and write a table by hand....

Regarding the other questions, I'll first have to have a look to the thread from March and the other suggestions.

Thanks,
Martina

mlennert · June 21, 2008, 12:04pm

On 20/06/08 15:33, Martina Schaefer wrote:

What I want to do is the following:
I have satellite photos from a glacier and we want to get an idea of the distribution of the crevasses over this glacier.
In a first step we want to assign every 100x100m grid cell of this glacier a value reflecting if there are very little, some, a lot of crevasses. As we don't know for the moment how to do this automatically (perhaps it would be possible to use the color in a next step), the idea was to define this value by "eyeballing". Any ideas how to do this a better way are welcome, too, but I think it's not easy.

Now the questions is "only" how to handle the grid and my values in Grass ...

You could use v.mkgrid to create a vector grid, add a column for number of crevasses with v.db.addcol, then v.digit this grid with background command 'd.rast YourSatelliteImage' and use the "Show attributes" function of v.digit to add the relevant number to the new column for each cell in the grid (you have to click on the centroid of the cell to see the attributes).

Moritz

Martina_Schaefer · June 25, 2008, 1:15pm

Hi Moritz,

what you propose sounds good for me ... but it didn't work

What I want to do is the following:
I have satellite photos from a glacier and we want to get an idea of the distribution of the crevasses over this glacier.
In a first step we want to assign every 100x100m grid cell of this glacier a value reflecting if there are very little, some, a lot of crevasses. As we don't know for the moment how to do this automatically (perhaps it would be possible to use the color in a next step), the idea was to define this value by "eyeballing". Any ideas how to do this a better way are welcome, too, but I think it's not easy.

Now the questions is "only" how to handle the grid and my values in Grass ...

You could use v.mkgrid to create a vector grid,

I created such a grid. This grid does not appear in the list ¨Select item" when I do "add a raster layer" in the GIS Manager. I hope this is normal?
As Rows and Columns I put the rows an columns of my x-y gridding, I hope this is correct?

add a column for number of crevasses with v.db.addcol,

in v.db.addcol my new grid appears in the list of items for "vector map for which to edit attribut table), I put layer=1 and name=D
then I get

"DBMI-DBF driver error:
SQL parser error in statement:
ALTER TABLE density ADD COLUMN D
Error in db_execute_immediate()
Error while executing: "ALTER TABLE density ADD COLUMN D
"
ERROR: Cannot continue (problem adding column)."

Has anybody an idea what this could mean?

Thanks,
Martina

Martina_Schaefer · June 25, 2008, 1:49pm

Hi again,

What I also want to do, is measure distances and angles between
objects on my images. For the moment I don't really see where
to search for this kind of tools. I thought that it would be possible
to draw lines and measure them and angles between them?

m.cogo might help. or 'v.distance upload=to_along,to_angle', d.measure, v.to.db...

In general beware to double check results/meaning if using a lat/lon location (see d.rhumbline and d.geodesic), or map projection with possibly heavy distortion at far distances.

If you build the Python SWIG interface there is also a m.distance module to look at.

see also this thread+code patch from March:
"How to find out an angle between to points on the map"
http://thread.gmane.org/gmane.comp.gis.grass.user/22817/
http://thread.gmane.org/gmane.comp.gis.grass.devel/25396/

d.measure works fine finally for measuring distances (even if I don't get to work the command line d.measure but only the one in menu of the Map Display).
But I still don't have an idea about how to measure the angle between two of this drawn lines.
I went through the thread indicated about, but I didn't understand anything :-(.

Thanks for your help,
Martina

mlennert · June 25, 2008, 2:33pm

On 25/06/08 15:15, Martina Schaefer wrote:

Hi Moritz,

what you propose sounds good for me ... but it didn't work

What I want to do is the following:
I have satellite photos from a glacier and we want to get an idea of the distribution of the crevasses over this glacier.
In a first step we want to assign every 100x100m grid cell of this glacier a value reflecting if there are very little, some, a lot of crevasses. As we don't know for the moment how to do this automatically (perhaps it would be possible to use the color in a next step), the idea was to define this value by "eyeballing". Any ideas how to do this a better way are welcome, too, but I think it's not easy.

Now the questions is "only" how to handle the grid and my values in Grass ...

You could use v.mkgrid to create a vector grid,

^^^^^^

I created such a grid. This grid does not appear in the list ¨Select item" when I do "add a raster layer" in the GIS Manager. I hope this is normal?

Yes, see hint just above

As Rows and Columns I put the rows an columns of my x-y gridding, I hope this is correct?

What do you mean by "my x-y gridding".

As you want to have grid cells of a given size, you will have to calculate your number of grids.

You have two options for creating your grid:

1) using position=region
- Set your region extents to multiples of your grid resolution
(i.e. 100m in your case) - the origin of your grid will be the lower left corner defined by your region
- Calculate your numbers of columns and rows using the region extension info (g.region -p), i.e. rows=North-South/Resolution, cols=East-West/Resolution

2) using position=coor (idependent of region settings)
  - Choose a point at the South-West corner just outside of your
image as starting point. Depending on your needs you might want to round these coordinates to your resolution.
  - Use r.info on your image to identify its extents
  - Calculate your numbers of columns and rows using the image extents, i.e. rows=North-South/Resolution, cols=East-West/Resolution

add a column for number of crevasses with v.db.addcol,

in v.db.addcol my new grid appears in the list of items for "vector map for which to edit attribut table), I put layer=1 and name=D
then I get

"DBMI-DBF driver error:
SQL parser error in statement:
ALTER TABLE density ADD COLUMN D
Error in db_execute_immediate()
Error while executing: "ALTER TABLE density ADD COLUMN D
"
ERROR: Cannot continue (problem adding column)."

Has anybody an idea what this could mean?

Yes, you need to indicate the type of column as well, so the columns parameter should be 'D int' if D is supposed to hold the number of crevasses which I believe would be integers.

The general idea of my proposal was to create a _vector_ grid with the relevant column in its attribute table and then use the digitization tool v.digit to display this grid on top of your image and edit the attribute values after visually identifying the number of crevasses in your grid cell. You would then end up with a vector grid map with the crevasse count in each cell. If necessary, you can transform this into raster with v.to.rast.

Another option to replace v.digit would be to use the command line and "old-style" monitors to display your two layers and then use d.what.vect -e to edit the attributes. Something like this (in a terminal - after having added your count column to the vector grid's attribute table):

d.mon x0
d.rast YourImage
d.vect YourGrid type=area fcol=none
d.what.vect -e map=YourGrid (-> click on the map to see a form pop up)

Moritz

Hamish · June 25, 2008, 4:41pm

Martina:

>> What I also want to do, is measure distances and
>> angles between objects on my images. For the moment I don't
>> really see where to search for this kind of tools. I thought that
>> it would be possible to draw lines and measure them and angles
>> between them?

...

d.measure works fine finally for measuring distances (even
if I don't get to work the command line d.measure but only
the one in menu of the Map Display).
But I still don't have an idea about how to measure the
angle between two of this drawn lines.

--this won't work correctly for lat/lon, only projected locations--
maybe the easies way is to use d.where (use middle mouse button to draw a line)

then take the two east and north coordinates:

# pythag hypot: a^2 + b^2 = c^2
distance = sqrt( (e2-e1)^2 + (n2-n1)^2 )

# angle: tan(theta) = opposite side of triange / adjacent side
theta = atan2(n2-n1, e2-e1)

$radians2deg:
theta_degrees = theta * 180/pi

# convert angle counterclockwise-from-east to clockwise-from-north:
heading = 90 - theta_degrees
if(heading < 0) then heading = heading + 360

you could write a little script to process all that from two mouse clicks.
d.where | awk ...
or
COORDS=`d.where`
.....

but that is left as an exercise for another day.

Hamish

Martina_Schaefer · June 25, 2008, 5:14pm

Hi,

thanks again.

Hamish wrote:

Martina:

What I also want to do, is measure distances and
angles between objects on my images. For the moment I don't
really see where to search for this kind of tools. I thought that
it would be possible to draw lines and measure them and angles
between them?

...

d.measure works fine finally for measuring distances (even
if I don't get to work the command line d.measure but only
the one in menu of the Map Display).
But I still don't have an idea about how to measure the
angle between two of this drawn lines.

--this won't work correctly for lat/lon, only projected locations--

Yes, I need the lat/lon files only to print some maps for somebody and I'm doing (or at least trying) to do my analysis on the UTM files.

maybe the easies way is to use d.where (use middle mouse button to draw a line)
then take the two east and north coordinates:
# pythag hypot: a^2 + b^2 = c^2
distance = sqrt( (e2-e1)^2 + (n2-n1)^2 )
# angle: tan(theta) = opposite side of triange / adjacent side
theta = atan2(n2-n1, e2-e1)
$radians2deg:
theta_degrees = theta * 180/pi
# convert angle counterclockwise-from-east to clockwise-from-north:
heading = 90 - theta_degrees if(heading < 0) then heading = heading + 360
you could write a little script to process all that from two mouse clicks.
d.where | awk ...
or
COORDS=`d.where`
.....

ok,... so there is nothing implemented
I'll see, either I'll write this script or I'll just print a map and measure the angles by hand, should be faster and good enough for what we need!

Martina

Hamish · June 26, 2008, 2:57am

Martina

Yes, I need the lat/lon files only to print some maps for somebody
and I'm doing (or at least trying) to do my analysis on the
UTM files.

note you can overlay a lat/lon grid in the UTM projection with
d.grid -g or -w, or ps.map's geogrid instruction.

Hamish

Martina_Schaefer · July 3, 2008, 4:17pm

Hi Hamish,

Yes, I need the lat/lon files only to print some maps for somebody
and I'm doing (or at least trying) to do my analysis on the
UTM files.

note you can overlay a lat/lon grid in the UTM projection with d.grid -g or -w, or ps.map's geogrid instruction.

Thanks for this hint, it's indeed what I wanted to do ... in the meanwhile I managed to work on another computer with ArcGIS and I could print all my maps.
With Grass I still have a lot of troubles...
Regarding the d.grid, I can't chose for the moment different gridspaces in lot and lat what I would like to.
And overall, I have not yet really understood how to deal with the graphicsmonitor.
Normally I do my plots just by clicking with the GIS Manager and then everything appears in 'Map Display 1'. But when I start a graphics monitor with d.mon x0 I don't manage to get everything in this window, for ex. my sat. data here is a RGB layer data, with d.rast I get each color as black and white in the x0 window, but not my real map....

Perhaps I didn't look at the right place, but up to now I didn't found a documentation on Grass that helped me with this and overall to find commands adapted to what I want to do...

Thanks for your help anyway

Martina