[GRASSLIST:9180] [Re: GRASS 6 in a nutshell tutorial]

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected California layout shapefile (form <www.census.gov>) and then importing it via v.proj in the location with the digital elevation map. Apparently the official shapefile of California layout has a wrong latitude.

Next step is to connect the output table of the population model (population densities at weather stations) to the map I have just created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth to manually put the file in that format or do you think is better to convert it in a database to connect with a specific command?

Thank you,

Luigi

-------- Original Message --------
Subject: [GRASSLIST:9087] Re: GRASS 6 in a nutshell tutorial
Date: Wed, 16 Nov 2005 00:03:56 +0100
From: lponti@infinito.it
To: Ian MacMillan <ian_macmillan@umail.ucsb.edu>
CC: Hamish <hamish_nospam@yahoo.com>, GRASSLIST@baylor.edu
References: <43667143.6030900@infinito.it> <43689B17.9040302@itc.it> <436AA43E.2030709@infinito.it> <2a710737c668a9cf27237f3236b55a2f@umail.ucsb.edu> <43724994.2070505@infinito.it> <20051110181214.21c19d6a.hamish_nospam@yahoo.com> <43790A77.2050108@infinito.it> <29ded48ac16bbfc3492e5eff023d357a@umail.ucsb.edu>

Dear Ian,

the original metadata of my state outline shapefile is:

PROJCS["unnamed",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG",7019]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6269]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat","NORTH"],AXIS["Long","EAST"],AUTHORITY["EPSG",4269]],PROJECTION["Albers"],PARAMETER["standard_parallel_1",34],PARAMETER["standard_parallel_2",40.5],PARAMETER["latitude_of_center",0],PARAMETER["Central_Meridian",-120],PARAMETER["False_Easting",0],PARAMETER["False_Northing",-4000000],UNIT["m",1]]

Based on this, following direction in the GRASS/cygwin installation page
http://geni.ath.cx/grass.html#toc15

I have created a dummy vector for reprojection of the raster DEM of western USA - e10g tile. Then I created a new location from EPSG code to match the vector .shp state outline file, and I have imported both the dummy raster DEM via "v.proj" and the vector state outline via "v.in.ogr".

Using "gdalinfo" on the original raster file (e10g), originated the following output:

gdalinfo c:/cygwin/home/andy/DEM_data/e10g
Driver: EHdr/ESRI .hdr Labelled
Size is 10800, 6000
Coordinate System is `'
Origin = (-180.000000,50.000000)
Pixel Size = (0.00833333,-0.00833333)
Corner Coordinates:
Upper Left (-180.0000000, 50.0000000)
Lower Left (-180.0000000, 0.0000000)
Upper Right ( -90.0000000, 50.0000000)
Lower Right ( -90.0000000, 0.0000000)
Center (-135.0000000, 25.0000000)
Band 1 Block=10800x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-500

Do you see anything wrong?

In addition, I tried to import v.in.e00 (after adding the required .exe file in the bin folder of GRASS) but I have not been successful yet.

Regards and thanks,

Luigi

On Mon, 14 Nov 2005 18:31:43 -0800
Ian MacMillan <ian_macmillan@umail.ucsb.edu> wrote:

Luigi, if you look at r.info and v.info, you will only find the projection information of your location, not necessarily the projection/datum/etc. information of the original data (I believe). Do you have the original metadata? If not, try running gdalinfo on your original data files to determine their original projection info.

man gdalinfo
will give you the correct syntax.

G'luck
Ian

On Nov 14, 2005, at 2:06 PM, Luigi Ponti wrote:

Thanks Amish,

looking at v.info and r.info, it seems the two maps have the same projection (see below). However, when looking at the output "d.where -w" for a sample map point (the San Francisco Bay), the longitude is very different in the two maps:

                                                  WGS84 Co-ordinates
            EAST: NORTH: LON:
            -191114.6346037 -24240.66895625 143:03:58.588289W 29:24:27.1
(San Francisco Bay on the California State vector layer - WRONG)

1839139.30311877 663981.00484798 122:19:50.783112W 37:57:48.1
(San Francisco Bay on the Elevation Model raster layer - CORRECT)

Do you think this could be due to an error in the import process or should I just look for another vector map?

Thanks again,

Luigi

v.info map=state layer=1
+---------------------------------------------------------------------- ------+
| Layer: state Organization:
                    |
| Mapset: luigi Source Date:
                     |
| Location: Cal_nad83-zone3 Name of creator: |
| Database: /home/andy
                                                      |
| Title:
                                                          
         |
| Map Scale: 1:1
                                                           |
| Map format: native
                                                        |
|---------------------------------------------------------------------- ------|
| Type of Map: Vector (level: 2)
                                         |
|
                                                          
                |
| Number of points: 0 Number of areas: 51 |
| Number of lines: 0 Number of islands: 48 |
| Number of boundaries: 59 Number of faces: 0 |
| Number of centroids: 51 Number of kernels: 0 |
|
                                                          
                |
| Map is 3D: 0
                                               |
| Number of dblinks: 1
                                               |
|
                                                          
                |
| Projection: Lambert Conformal Conic (zone 0)
                            |
| N: 922206.900 S: -701781.905
                                |
| E: 1043412.917 W: -374425.808
                              |
| B: 0.000 T: 0.000
                                      |
|
                                                          
                |
| Digitize threshold: 0.00000
                                             |
| Comments:
                                                               |
|
                                                          
                |
+---------------------------------------------------------------------- ------+

r.info map=e10g@PERMANENT
+---------------------------------------------------------------------- ------+
| Layer: e10g Date: Tue Nov 8 17:25:24 2005 |
| Mapset: PERMANENT Login of Creator: andy |
| Location: Cal_nad83-zone3
                                                 |
| DataBase: /home/andy
                                                      |
| Title: ( e10g )
                                                       |
|---------------------------------------------------------------------- ------|
|
                                                          
                |
| Type of Map: raster Number of Categories: 5443 |
| Data Type: CELL
                                                      |
| Rows: 7192
                                                      |
| Columns: 11628
                                                     |
| Total Cells: 83628576
                                                  |
| Projection: Lambert Conformal Conic (zone 0)
                       |
| N: 3353466 S: -3838238 Res: 999.95884316 |
| E: 6073284 W: -5555183 Res: 1000.04016168 |
| Range of data: min = -84 max = 5443
                                 |
|
                                                          
                |
| Data Source:
                                                            |
|
                                                          
                |
|
                                                          
                |
|
                                                          
                |
| Data Description:
                                                       |
| generated by r.proj
                                                    |
|
                                                          
                |
|
                                                          
                |
+---------------------------------------------------------------------- ------+

Hamish wrote:

I have imported both the vector and the raster with minimum trouble, since the first has a .prj file which GRASS recognized, and the second

had a detailed import tutorial in the cygwin-GRASS install page. However, despite the fact the the two layers supposedly have the same NAD83 projection info (I did not enter anything manually), when I try
to display them together on a GRASS monitor, the state layout is in
the middle of the ocean compared to the elevation raster.

Is there a way to overlap the two layers or to chek for a coordinate
error?

r.info and v.info will give you the map extents, you should then be able
to figure out which is correct.

Also try "d.where -w" to give you a lat/lon value for different parts of your map which you can compare to the map on the wall.

Hamish

>
What happens if a big asteroid hits Earth? Judging from realistic simulations involving a sledge hammer and a common laboratory frog, we can assume it will be pretty bad.
- Dave Barry

-------------------------------------------------------------
This message has been scanned by Postini anti-virus software.

INFINITO ADSLFLAT 4 MEGA: SOLO 27,90 EURO AL MESE IVA INCLUSA IP STATICO, BANDA GARANTITA 256Kbps ANTIVIRUS E FIREWALL INCLUSI NEL PREZZO

http://adsl.infinito.it

Hello Luigi,

On Mon, 21 Nov 2005 17:20:05 -0800 Luigi Ponti <lponti@infinito.it>
wrote:

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected
California layout shapefile (form <www.census.gov>) and then
importing it via v.proj in the location with the digital elevation
map. Apparently the official shapefile of California layout has a
wrong latitude.

Next step is to connect the output table of the population model
(population densities at weather stations) to the map I have just
created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with
v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

you could use the fs="," for selecting a , as a selector.

have a look at the manual-page of v.in.ascii.

g.manual v.in.ascii or

v.in.ascii -h

Best
  Stephan

--
GDF Hannover - Solutions for spatial data analysis and remote sensing
Hannover Office - Mengendamm 16d - D-30177 Hannover
Internet: www.gdf-hannover.de - Email: holl@gdf-hannover.de
Phone : ++49-(0)511.39088507 - Fax: ++49-(0)511.39088508

Thanks Stephan,

I was able to import the ASCII text file to vector points and interpolate between them.
Also, based on a very useful but unrelated discussion, I obtained a very pretty shaded relief over which I draped the interpolated raster.

Now I have started clipping the region that are not of interest to the analysis. To clip out what is not California, I have transformed a vector layout of the state to a raster via v.to.rast, and used it as a mask in v.surf.rst.

How could I clip regions that have an elevation greater than a particular value based on the elevation model raster?

Thanks again,

Luigi

Stephan Holl wrote:

Hello Luigi,

On Mon, 21 Nov 2005 17:20:05 -0800 Luigi Ponti <lponti@infinito.it>
wrote:

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected California layout shapefile (form <www.census.gov>) and then
importing it via v.proj in the location with the digital elevation
map. Apparently the official shapefile of California layout has a
wrong latitude.

Next step is to connect the output table of the population model (population densities at weather stations) to the map I have just created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with
v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

you could use the fs="," for selecting a , as a selector.

have a look at the manual-page of v.in.ascii.

g.manual v.in.ascii or

v.in.ascii -h

Best
Stephan

Luigi, try something like:
r.mapcalc newmap= 'if (oldmap>1000,oldmap,null())'

-Ian

On Nov 28, 2005, at 1:48 PM, Luigi Ponti wrote:

Thanks Stephan,

I was able to import the ASCII text file to vector points and interpolate between them.
Also, based on a very useful but unrelated discussion, I obtained a very pretty shaded relief over which I draped the interpolated raster.

Now I have started clipping the region that are not of interest to the analysis. To clip out what is not California, I have transformed a vector layout of the state to a raster via v.to.rast, and used it as a mask in v.surf.rst.

How could I clip regions that have an elevation greater than a particular value based on the elevation model raster?

Thanks again,

Luigi

Stephan Holl wrote:

Hello Luigi,
On Mon, 21 Nov 2005 17:20:05 -0800 Luigi Ponti <lponti@infinito.it>
wrote:

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected California layout shapefile (form <www.census.gov>) and then
importing it via v.proj in the location with the digital elevation
map. Apparently the official shapefile of California layout has a
wrong latitude.

Next step is to connect the output table of the population model (population densities at weather stations) to the map I have just created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with
v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

you could use the fs="," for selecting a , as a selector.

have a look at the manual-page of v.in.ascii.

g.manual v.in.ascii or

v.in.ascii -h

Best
  Stephan

>
What happens if a big asteroid hits Earth? Judging from realistic simulations involving a sledge hammer and a common laboratory frog, we can assume it will be pretty bad.
  - Dave Barry

-------------------------------------------------------------
This message has been scanned by Postini anti-virus software.

Thanks Ian,

I used it and it worked very well. I also used only the r.mapcalc to clip instead of MASK, so as to clip just one layer and not the entire display.

I would like to ask if somebody knows the reason why when saving a display to a .png file with d.out.png, the legend produced via d.legend does not show up in the saved .png image.

As a side note, I would like to thank again the mailing: my mapping task is pretty much accomplished.

Luigi

Ian MacMillan wrote:

Luigi, try something like:
r.mapcalc newmap= 'if (oldmap>1000,oldmap,null())'

-Ian

On Nov 28, 2005, at 1:48 PM, Luigi Ponti wrote:

Thanks Stephan,

I was able to import the ASCII text file to vector points and interpolate between them.
Also, based on a very useful but unrelated discussion, I obtained a very pretty shaded relief over which I draped the interpolated raster.

Now I have started clipping the region that are not of interest to the analysis. To clip out what is not California, I have transformed a vector layout of the state to a raster via v.to.rast, and used it as a mask in v.surf.rst.

How could I clip regions that have an elevation greater than a particular value based on the elevation model raster?

Thanks again,

Luigi

Stephan Holl wrote:

Hello Luigi,
On Mon, 21 Nov 2005 17:20:05 -0800 Luigi Ponti <lponti@infinito.it>
wrote:

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected California layout shapefile (form <www.census.gov>) and then
importing it via v.proj in the location with the digital elevation
map. Apparently the official shapefile of California layout has a
wrong latitude.

Next step is to connect the output table of the population model (population densities at weather stations) to the map I have just created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with
v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

you could use the fs="," for selecting a , as a selector.

have a look at the manual-page of v.in.ascii.

g.manual v.in.ascii or

v.in.ascii -h

Best
    Stephan

>
What happens if a big asteroid hits Earth? Judging from realistic simulations involving a sledge hammer and a common laboratory frog, we can assume it will be pretty bad.
- Dave Barry

-------------------------------------------------------------
This message has been scanned by Postini anti-virus software.

Luigi, I think that I have seen on this list that you need to place your legend without using the mouse. Maybe you could place it on the screen using screen coordinates or map coordinates instead?

G'luck
-Ian

On Nov 29, 2005, at 1:11 PM, Luigi Ponti wrote:

Thanks Ian,

I used it and it worked very well. I also used only the r.mapcalc to clip instead of MASK, so as to clip just one layer and not the entire display.

I would like to ask if somebody knows the reason why when saving a display to a .png file with d.out.png, the legend produced via d.legend does not show up in the saved .png image.

As a side note, I would like to thank again the mailing: my mapping task is pretty much accomplished.

Luigi

Ian MacMillan wrote:

Luigi, try something like:
r.mapcalc newmap= 'if (oldmap>1000,oldmap,null())'

-Ian

On Nov 28, 2005, at 1:48 PM, Luigi Ponti wrote:

Thanks Stephan,

I was able to import the ASCII text file to vector points and interpolate between them.
Also, based on a very useful but unrelated discussion, I obtained a very pretty shaded relief over which I draped the interpolated raster.

Now I have started clipping the region that are not of interest to the analysis. To clip out what is not California, I have transformed a vector layout of the state to a raster via v.to.rast, and used it as a mask in v.surf.rst.

How could I clip regions that have an elevation greater than a particular value based on the elevation model raster?

Thanks again,

Luigi

Stephan Holl wrote:

Hello Luigi,
On Mon, 21 Nov 2005 17:20:05 -0800 Luigi Ponti <lponti@infinito.it>
wrote:

Thanks Ian and Hamish,

I was eventually able to fix the problem by downloading a unprojected California layout shapefile (form <www.census.gov>) and then
importing it via v.proj in the location with the digital elevation
map. Apparently the official shapefile of California layout has a
wrong latitude.

Next step is to connect the output table of the population model (population densities at weather stations) to the map I have just created. This would hopefully result in an interpolated raster.

Hamish suggested the following:

Load in site data with v.in.ascii & interpolate between sites with
v.surf.rst.

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

you could use the fs="," for selecting a , as a selector.

have a look at the manual-page of v.in.ascii.

g.manual v.in.ascii or

v.in.ascii -h

Best
    Stephan

>
What happens if a big asteroid hits Earth? Judging from realistic simulations involving a sledge hammer and a common laboratory frog, we can assume it will be pretty bad.
- Dave Barry

-------------------------------------------------------------
This message has been scanned by Postini anti-virus software.

>
What happens if a big asteroid hits Earth? Judging from realistic simulations involving a sledge hammer and a common laboratory frog, we can assume it will be pretty bad.
  - Dave Barry

-------------------------------------------------------------
This message has been scanned by Postini anti-virus software.

Yes, that was it! Placing the legend on the screen using screen coordinates worked, as the legend displayed correctly in the .png file.

I am now trying to import a shapefile (with .dbf and .shx dependencies in the same folder). It contains a Reference EvapoTranspiration Zones Map for California that I have obtained by contacting California Irrigation Management Information System (CIMIS) <http://wwwcimis.water.ca.gov/&gt;\. I am using v.in.ogr to import the shapefile to a latlong location by overriding the projection because I don't know anything about it.
This is the shell output I get:

[output starts]

GRASS 6.1.cvs (latlong):~ > v.in.ogr dsn=/home/Administrator/GISBASE/CIMIS/ outp
ut=eto3 layer=eto_utm10_nad27 min_area=0.0001 snap=-1 -o
Over-riding projection check.
Proceeding with import...
Layer: eto_utm10_nad27
WARNING: Area size 0.0e+00, area not imported.

... omissis ... (about 40 lines of the same warning: "WARNING: Area size 0.0e+00, area not imported.")

WARNING: Area size 0.0e+00, area not imported.
-----------------------------------------------------
Building topology ...
344 primitives registered Building areas: 100%
342 areas built 339 isles built
Attaching islands: 100%
Attaching centroids: 100%
Topology was built.
Number of nodes : 341
Number of primitives: 344
Number of points : 0
Number of lines : 0
Number of boundaries: 344
Number of centroids : 0
Number of areas : 342
Number of isles : 339
Number of incorrect boundaries : 4
Number of areas without centroid : 342
-----------------------------------------------------
WARNING: Cleaning polygons, result is not guaranteed!
Building topology ...
Topology was built.
Number of nodes : 341
Number of primitives: 344
Number of points : 0
Number of lines : 0
Number of boundaries: 344
Number of centroids : 0
Number of areas : -
Number of isles : -
-----------------------------------------------------
Break polygons:
Registering points ... 3372
All points (vertices): 4245
Registered points (unique coordinates): 3372
Points marked for break: 450
Breaks: 512
-----------------------------------------------------
Remove duplicates:
Duplicates: 133
-----------------------------------------------------
Break boundaries:
Intersections: 0 -----------------------------------------------------
Remove duplicates:
Duplicates: 0
-----------------------------------------------------
Clean boundaries at nodes:
Modifications: 0
-----------------------------------------------------
Change dangles to lines:
Removed dangles: 0 removed lines: 0
-----------------------------------------------------
Remove bridges:
Removed bridges: 0 removed lines: 0
-----------------------------------------------------
Building topology ...
Building areas: 100%
353 areas built 216 isles built
Attaching islands: 100%
Topology was built.
Number of nodes : 451
Number of primitives: 857
Number of points : 0
Number of lines : 0
Number of boundaries: 857
Number of centroids : 0
Number of areas : 353
Number of isles : 216
Number of areas without centroid : 353
Layer: eto_utm10_nad27
-----------------------------------------------------
Building topology ...
931 primitives registered Building areas: 100%
353 areas built 216 isles built
Attaching islands: 100%
Attaching centroids: 100%
Topology was built.
Number of nodes : 794
Number of primitives: 931
Number of points : 0
Number of lines : 0
Number of boundaries: 587
Number of centroids : 344
Number of areas : 353
Number of isles : 216
Number of areas without centroid : 9
-----------------------------------------------------
392 input polygons
total area: 2.689311e+16 (353 areas)
overlapping area: 0.000000e+00 (0 areas)
area without category: 1.036357e+15 (9 areas)

[output ends]

After importing the vector, it is not possible to display it ("The bounding box of the map outside current region, nothing displayed.") or even to g.remove it. However, when I try to open it with QGIS, it displays with no problem and gives the following "Layer Spatial Reference System": +proj=longlat +ellps=WGS84 +no_defs.

I also tried to import the vector to a location with "utm10_nad27"-like projection, as the name of the file would suggest. In this case, GRASS complains about a coor file which is too big.

If anybody has suggestions, please let me know.

Thanks,

Luigi

Ian MacMillan wrote:

Luigi, I think that I have seen on this list that you need to place your legend without using the mouse. Maybe you could place it on the screen using screen coordinates or map coordinates instead?

G'luck
-Ian

On Nov 29, 2005, at 1:11 PM, Luigi Ponti wrote:

Thanks Ian,

I used it and it worked very well. I also used only the r.mapcalc to clip instead of MASK, so as to clip just one layer and not the entire display.

I would like to ask if somebody knows the reason why when saving a display to a .png file with d.out.png, the legend produced via d.legend does not show up in the saved .png image.

As a side note, I would like to thank again the mailing: my mapping task is pretty much accomplished.

Luigi

Ian MacMillan wrote:

After importing the vector, it is not possible to display it ("The
bounding box of the map outside current region, nothing displayed.")

Use g.region to zoom to the extents of the imported map first,

g.region vect=map
d.vect map
v.info map

Hamish

According to the manual v.in.ascii uses the | separator. Is it worth
to manually put the file in that format or do you think is better to
convert it in a database to connect with a specific command?

You can use the v.in.ascii fs= option to set the separator to whatever
you like. e.g. "fs=," or "fs=tab"

Hamish