Re: [GRASS-user] v.what.strds error creating column with @ name
Hi Markus,
>> I'm trying to work in multiple mapsets and with multiple tables. I
>> find this hard to do when they're in different databases, so when I
>> create a new mapset I set the DB connection to PERMANENT:
>>
>> db.connect database=${LOCATION_NAME}/PERMANENT/sqlite/sqlite.db
>
> I believe that the quotes are missing - so the variable is immediately
> interpreted.
I don't think it matters? With quotes it will resolve to real paths sometime
later, but LOCATION_NAME never changes, so when/where the variable expansion
occurs shouldn't change anything. If $MAPSET were in quotes, then it would
create a new DB per mapset, which I do not want, hence hard-coding PERMANENT.
I can't read the URL you sent - I have email but not web browser for another
week. But on my local computer grass-sqlite.html manual page says
> [...] the file storage location can be freely chosen.
Suggesting I could hard-code this anywhere I want, even outside of the GRASS
folder structure.
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.what.strds error creating column with @ name
Hi Veronica,
On 2024-02-03 at 02:52 +13, Veronica Andreo
wrote...
> Were you able to solve the problem or find its cause? Could you create
> a reproducible example with the NC dataset? It seems really strange,
> if the mapset is named vel why would it add the date to it too?
I'm guessing the problem and cause is me, not GRASS. Before I try to replicate
in the NC data set, I hope you can help answer a question.
I'm trying to work in multiple mapsets and with multiple tables. I find this
hard to do when they're in different databases, so when I create a new mapset I
set the DB connection to PERMANENT:
db.connect database=${LOCATION_NAME}/PERMANENT/sqlite/sqlite.db
I assume this may be part of the problem. Is this common practice? Or is there
some other way to work with `db.execute` and multiple databases?
Thanks,
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] v.what.strds error creating column with @ name
Hi List, I'm trying to sample an STRDS at some vector points. My points are generated (from lines in another mapset) with: v.to.points input=ex1@gates type=line output=points dmax=100 My STRDS is generated in the current mapset (named vel), and then when I run v.what.strds input=points strds=vx,vy output=test I get the following error: DBMI-SQLite driver error: Error in sqlite3_prepare(): near "@vel_2007_09_07": syntax error ERROR: Error while executing: 'ALTER TABLE test_points_1 ADD COLUMN vy_2007_09_07_2008_04_23@vel_2007_09_07 DOUBLE PRECISION' ERROR: Unable to add column . ERROR: Unable to add column vy_2007_09_07_2008_04_23@vel_2007_09_07 DOUBLE PRECISION to vector map '@' is not a valid character name in a DB column. Am I doing something wrong that v.what.strds is trying to generate this column name? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to replace some vector & DB features, keep attributes
Hi Maris, On 2024-01-15 at 02:10 -08, Maris Nartiss wrote... > v.edit tool=catadd/catdel is the thing you are looking for. Delete old > geometry, delete category of new geometry (if there is one), set old > cat value to the new geometry. Thank you for the suggestion. I don't need cat to be the same. I'm deleting ~100 cats (old vector lines) and replacing with only one new one anyway. More importantly, in my simplified example I neglected to mention that although I want to keep some attributes, I don't want to keep all. There are several old attributes that are no longer relevant for this one manually created line. I was going to create the new line (as below) and then in a for loop: for a in attr1 attr2 attr3 attr4; do # attrs I want to keep or copy # db.select statement to get attr1 value from $old_cat # db.execute sql='UPDATE test set a = $att where cat=$new_cat done Alternatively, if I follow your advice, then I'd need to: for a in attr1 attr2 attr3 attr4; do # attrs I want to delete # db.execute sql='UPDATE test set $attr = NULL where cat=$new_cat done So... is this looping method correct, or is there some SQL statement that will copy attr1,attr2,attrn from row x (defined by an ID or cat) to row y? Thanks, -k. >> I'd like to replace multiple vector lines with ID == 42 with one new >> one, and keep the attributes (date, year, name) from one of the >> replaced/removed lines. >> >> I'm able to create my one new line with: >> >> >> I'm not sure how to copy over the relevant attributes from one of the >> lines (first, last, doesn't matter) that has the relevant ID. >> >> Once I do that, I think I know how to delete the lines I don't want, >> but I'm not sure I'm doing it correct. It seems I have to run two >> commands to delete the lines from both the vectors (when displaying) >> and the database (when querying). I'm doing this: >> >> v.edit map=test tool=delete where="ID == 42" db.execute sql="DELETE >> FROM test WHERE ID == 42" >> >> Can anyone help with the middle step? >> >> Thanks, >> ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] How to replace some vector & DB features, keep attributes
Hi, I'd like to replace multiple vector lines with ID == 42 with one new one, and keep the attributes (date, year, name) from one of the replaced/removed lines. I'm able to create my one new line with: echo "L 2 1 -273157 -994455 -255458 -989423 1 10" | v.edit -n tool=add map=test input=- I'm not sure how to copy over the relevant attributes from one of the lines (first, last, doesn't matter) that has the relevant ID. Once I do that, I think I know how to delete the lines I don't want, but I'm not sure I'm doing it correct. It seems I have to run two commands to delete the lines from both the vectors (when displaying) and the database (when querying). I'm doing this: v.edit map=test tool=delete where="ID == 42" db.execute sql="DELETE FROM test WHERE ID == 42" Can anyone help with the middle step? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Merging stream segments based on order
Please disregard last message. The 'merge' tool example uses 'cat', but realized I can select on any table column, and the tool takes care of the multiple unconnected (different) streams/basins. v.edit map=streams tool=merge where="hack == 1" v.edit map=streams tool=merge where="hack == 2" v.edit map=streams tool=merge where="hack == 3" and then calculate sinuosity with v.to.db map=streams option=sinuous columns=sinuosity Seems to work. Thank you, -k. On 2024-01-02 at 13:53 -07, Ken Mankoff wrote... > Hi GRASS list, > > I'd like to calculate stream sinuosity. I currently do this with: > > v.to.db map=streams option=sinuous columns=sinuosity > > But the 'streams' vector there is generated from > > r.stream.order stream_rast=streams direction=fdir elevation=head > accumulation=acc stream_vect=streams hack=hack shreve=shreve --o > > And is made up of many segments. I'm therefore calculating sinuosity > of each segment. I believe this underestimates sinuosity. For example, > a stream can have many short straight lines (sinuosity 1), each > connected at an angle (sinuosity >1). But the reported sinuosity is > only 1 for each of the straight segments. > > I'd like to merge all connected lines with say, hack order 1, and then > calculate sinuosity for one much longer segment. The stream vector has > many basins, so there are many hack=1 segments. The associated table > has 40k rows. If I limit to where hack==1, there are 6k rows. hack==2 > has 12k rows. > > Can anyone suggest how I merge? Looping over basin seems inefficient. > I see v.edit has a 'merge' tool, but it isn't clear to me how to find > the cats I want to merge? For each segment I have: > > cat|stream|next_stream|prev_str01|prev_str02|prev_str03|prev_str04|strahler|horton|shreve|hack|topo_dim|scheidegger|drwal_old|length|stright|sinosoid|cum_length|flow_accum|out_dist|source_elev|outlet_elev|elev_drop|out_drop|gradient|start_x|start_y|stop_x|stop_y|sinuosity|n_up > > I could add two new new columns that are the out_x,out_y of the > terminal outlet (direction = -1 from r.watershed), build a uniq ID > from that pair, set that to the category of all of the segments, and > then merge using that? > > Thanks for any suggestions, > > -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Merging stream segments based on order
Hi GRASS list, I'd like to calculate stream sinuosity. I currently do this with: v.to.db map=streams option=sinuous columns=sinuosity But the 'streams' vector there is generated from r.stream.order stream_rast=streams direction=fdir elevation=head accumulation=acc stream_vect=streams hack=hack shreve=shreve --o And is made up of many segments. I'm therefore calculating sinuosity of each segment. I believe this underestimates sinuosity. For example, a stream can have many short straight lines (sinuosity 1), each connected at an angle (sinuosity >1). But the reported sinuosity is only 1 for each of the straight segments. I'd like to merge all connected lines with say, hack order 1, and then calculate sinuosity for one much longer segment. The stream vector has many basins, so there are many hack=1 segments. The associated table has 40k rows. If I limit to where hack==1, there are 6k rows. hack==2 has 12k rows. Can anyone suggest how I merge? Looping over basin seems inefficient. I see v.edit has a 'merge' tool, but it isn't clear to me how to find the cats I want to merge? For each segment I have: cat|stream|next_stream|prev_str01|prev_str02|prev_str03|prev_str04|strahler|horton|shreve|hack|topo_dim|scheidegger|drwal_old|length|stright|sinosoid|cum_length|flow_accum|out_dist|source_elev|outlet_elev|elev_drop|out_drop|gradient|start_x|start_y|stop_x|stop_y|sinuosity|n_up I could add two new new columns that are the out_x,out_y of the terminal outlet (direction = -1 from r.watershed), build a uniq ID from that pair, set that to the category of all of the segments, and then merge using that? Thanks for any suggestions, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Convert vector to STVDS
Hi List, On 2023-12-01 at 16:30 -08, Ken Mankoff wrote... > I have a vector with 115284 lines. From db.select I see > > | cat | ID | Date | > | 1 | 278 | 1990-07-10 | > | 2 | 278 | 1992-07-30 | > | 3 | 278 | 1994-08-29 | > | 4 | 242 | 1998-06-05 | > | 5 | 255 | 2005-03-03 | > > That is, 115284 cats and 291 IDs. Using just one ID to start with (with 80 timestamps) I've created 80 vectors, each with their own layer and v.timestamp. I can then patch these together to one product with multiple layers. If I t.register all of these, I also see what I expect to see with g.gui.timeline. If I run "g.gui.animation stvds=t" it reports "Topology [...] is invalid". Is there a way to export STVDs and see the time component in QGIS or KMZ/Google Earth? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Convert vector to STVDS
Hello list, I have a vector with 115284 lines. From db.select I see | cat | ID | Date | | 1 | 278 | 1990-07-10 | | 2 | 278 | 1992-07-30 | | 3 | 278 | 1994-08-29 | | 4 | 242 | 1998-06-05 | | 5 | 255 | 2005-03-03 | | 6 | 255 | 1985-09-05 | | 7 | 255 | 1986-09-08 | | 8 | 278 | 1997-09-08 | That is, 115284 cats and 291 IDs. I'd like to work with this as an STVDS. I know how to go about this if it were all rasters, but am still (after years) wrapping my head around the GRASS vector format and power. Can someone advise how I would go about importing this into an STVDS? My first thought was to extract each line to a new vector. Obviously inefficient. Perhaps each ID to a new vector, and then each line (each unique Date) for each of those vectors is moved to a new layer? Thanks, -k. P.S. Unrelated but somewhat related. Does STR3DS let me work with 4D data: 3D (xyz) rasters that then have a time component added? Or am I still limited to 3D data? ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Probabilistic neighborhood analysis
What about using the fuzzy logic modules? -k. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Wed, Dec 14, 2022, 13:38 Bernardo Santos via grass-user < grass-user@lists.osgeo.org> wrote: > Hi, > > I am trying to produce scenarios of past land cover, before hydropower > reservoirs were built. To do so, I need to fill empty pixels from a raster > in the locations where the reservoirs are currently present, using as input > the actual land cover map. I tried doing that with r.neighbors (taking > method=mode) with neighborhoods of increasing size, to replace null pixels > with the most common land cover class in the neighborhood. I also tried > that with r.fill.stats which is basically the same thing. > However, the results gets very homogeneous, since the interpolated null > cells always get the value of the most common land cover class. > > Do anyway know of a method in GRASS to perform a "probabilistic" > neirighborhood analysis, where cells in a neighborhood are given weights > (possibly related to the distance to the central cell and to their > frequency) and these weights are used to stocastically sample a value to > fill the central cell? > If not in GRASS, does anyway know of such a method in a different > platform, i.e. R? > > Thanks! > Best > Bernardo > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user > ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] empty results from v.distance
Ah nevermind - I thought that all of the areas had values in the ID column, but that was incorrect, which is why the points vector was not fully populated. Thanks for helping me debug this. -k. On 2022-05-20 at 21:27 -07, Ken Mankoff wrote: > Hello GRASS list, > > I'm running v.distance because I'd like to assign the nearest area ID > to a points vector. Some of the points are outside any areas. > > I'm using this command: > > v.distance from=points to=areas upload=to_attr column=c to_column=ID > output=connections > > And I'd like to see a new column "C" in the points table that contains > the ID from the nearest area in areas. > > If I plot points, areas, and connections, it appears that every point > is connected to an area. But if I `db.select points|head`, many rows > in the ID column are empty. I would not expect any to be empty. Can > someone explain what I might be doing wrong? > > Thank you, > > -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] empty results from v.distance
Hello GRASS list, I'm running v.distance because I'd like to assign the nearest area ID to a points vector. Some of the points are outside any areas. I'm using this command: v.distance from=points to=areas upload=to_attr column=c to_column=ID output=connections And I'd like to see a new column "C" in the points table that contains the ID from the nearest area in areas. If I plot points, areas, and connections, it appears that every point is connected to an area. But if I `db.select points|head`, many rows in the ID column are empty. I would not expect any to be empty. Can someone explain what I might be doing wrong? Thank you, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] hexagon rasters
Hi Markus, Yes - that works but requires downstream hacks because of the different resolutions. Still, it is a workable solution I didn't think of. Thank you! -k. On 2022-02-09 at 12:52 -08, Markus Neteler wrote: > Hi Ken, > > On Tue, Feb 8, 2022 at 4:44 AM Ken Mankoff wrote: >> >> Hello List, >> >> I'm interested in working with hexagonal rasters. I know I can make >> hexagon vectors with "v.mkgrid -h", but is there any way to then work >> with these in raster space? I'd like to use r.walk and r.cost to >> estimate costs to move around the hex grid. > > I am sure you considered v.to.rast? > > # North Carolina sample dataset > g.region raster=elevation res=5000 -pa > > # create vector hexagons > v.mkgrid map=hexagons -h > v.info -c hexagons > > # convert to raster model, requires the resolution to be increased > (example: 5000m -> 250m) > g.region res=250 -p > v.to.rast input=hexagons output=hexagons use=cat > > # viz > d.rast hexagons > d.vect hexagons type=boundary > > Now you may turn them into cost surfaces. > > Best, > Markus ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] hexagon rasters
Hello List, I'm interested in working with hexagonal rasters. I know I can make hexagon vectors with "v.mkgrid -h", but is there any way to then work with these in raster space? I'd like to use r.walk and r.cost to estimate costs to move around the hex grid. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.fillnulls hang on one computer, not on another
Hi Markus, > I'm out of my depth here. It seems like a similar issue elsewhere may > point to SSL library issues or that this is a parent thread waiting on > another thread. From > https://meenakshi02.wordpress.com/2011/02/02/strace-hanging-at-futex/ > and 'ps -efL | grep filln', there is only one other PID. If I strace > that, all I see is > > read(3, It grows!! After a while I'm now seeing read(3, " 10%\10\10\10\10\10", 6410) = 10 read(3, " 20%\10\10\10\10\10", 6400) = 10 read(3, " 30%\10\10\10\10\10", 6390) = 10 read(3, " 40%\10\10\10\10\10", 6380) = 10 So maybe it'll finish... eventually? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.fillnulls hang on one computer, not on another
Hi Markus, On 2022-01-18 at 12:28 -08, Markus Neteler wrote: > On Tue, Jan 18, 2022 at 9:23 PM Ken Mankoff wrote: >> I'm trying to run 'r.fillnulls' and it hangs on one computer but not on >> another. >> > > and then run strace on it > > strace -p strace prints out things like this (at about 5 Hz): futex(0x55591e4fad24, FUTEX_WAIT_PRIVATE, 4832, NULL) = 0 futex(0x55591f6fcde4, FUTEX_WAKE_PRIVATE, 2147483647) = 16 futex(0x55591f6fcde4, FUTEX_WAIT_PRIVATE, 4928, NULL) = 0 futex(0x55591e4fad24, FUTEX_WAKE_PRIVATE, 2147483647) = 27 futex(0x55591f6fcde4, FUTEX_WAIT_PRIVATE, 4936, NULL) = -1 EAGAIN (Resource temporarily unavailable) futex(0x55591f6fcde4, FUTEX_WAIT_PRIVATE, 4952, NULL) = 0 futex(0x55591e4fad24, FUTEX_WAKE_PRIVATE, 2147483647) = 35 futex(0x55591f6fcde4, FUTEX_WAIT_PRIVATE, 4960, NULL) = 0 I'm out of my depth here. It seems like a similar issue elsewhere may point to SSL library issues or that this is a parent thread waiting on another thread. From https://meenakshi02.wordpress.com/2011/02/02/strace-hanging-at-futex/ and 'ps -efL | grep filln', there is only one other PID. If I strace that, all I see is read(3, And nothing else... -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.fillnulls hang on one computer, not on another
Hello, I'm trying to run 'r.fillnulls' and it hangs on one computer but not on another. On my laptop running Ubuntu 20.04 it runs fine in GRASS 7.8.6. It also runs in a docker container running GRASS 7.8.4. On our server it just hangs. This is running Ubuntu 18.04 and grass 7.8.2. It also hangs running *the same docker image* Ubuntu 20.04 GRASS 7.8.4 Do you have any suggestions how to debug this? Examining it with top sows it is churning away with 100 % CPU usage. It takes ~5 minutes on my laptop then completes. It never completes on the server (which is much more powerful - it should complete faster). Thank you, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Raster null values unchanged by r.null
I think r.null does not work on external (r.external) files. Make sure this is not the issue. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Wed, Oct 20, 2021, 07:02 Eric Patton via grass-user < grass-user@lists.osgeo.org> wrote: > I'm encountering some strange behaviour from r.null - when I use the > setnull parameter to assign a particular value to be null in a raster, the > raster areas just set to null are still visible and coloured according to > their previous values when I refresh the display in the gui map window. > Querying these values shows they haven't been set to null. > > Using version 7.8.5 on Linux Mint. Can anyone confirm? > > -- > Eric > > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user > ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Cannot read NetCDF with 'semi_major_axis' or 'inverse_flattening'
Ok - just to clear this up and shut down this thread. It turns out it is
EPSG:3411 but was reported as 3413. I can import (with -o) into a 3411
location, and then reproject.
Cheers,
-k.
On Sun, Jul 25, 2021 at 8:45 PM Ken Mankoff wrote:
> Hi Stefan,
>
> Thanks for the hint. I often use '-o'. In this case after discussion with
> the data provider it turns out that it was not EPSG:3413 but instead some
> custom projection, even though part of the NetCDF file reported EPSG:3413
> in the metadata.
>
> Unfortunately r.import does not work. I get the errors below. Given the
> 0..3..6..9... it appears to import something, but no raster exists when the
> command finishes.
>
> r.import input=NetCDF:file.nc:surface output=surface
> WARNING: Datum not recognised by GRASS and no parameters found
> WARNING: Projection of dataset does not appear to match current location.
>
> Location PROJ_INFO is:
> name: WGS 84 / NSIDC Sea Ice Polar Stereographic North
> datum: wgs84
> ellps: wgs84
> proj: stere
> lat_0: 90
> lat_ts: 70
> lon_0: -45
> x_0: 0
> y_0: 0
> no_defs: defined
> init: EPSG:3413
>
> Dataset PROJ_INFO is:
> name: unnamed
> a: 6378273
> es: 0.00669388300017
> proj: stere
> lat_0: 90
> lat_ts: 70
> lon_0: -45
> x_0: 0
> y_0: 0
> no_defs: defined
>
> ERROR: datum
>
> WARNING: Datum not recognised by GRASS and no parameters found
> WARNING: Datum not recognised by GRASS and no parameters found
> Importing raster map ...
>
> 0..3..6..9..12..15..18..21..24..27..30..33..36..39..42..45..48..51..54..57..60..63..66..69..72..75..78..81..84..87..90..93..96..99..100
> ERROR: Unable to open element file for
>
> ERROR: Unable to get reprojected map extent
>
>
> -k.
>
>
>
> On Sat, Jul 24, 2021 at 2:19 PM Stefan Blumentrath <
> stefan.blumentr...@nina.no> wrote:
>
>> Hei Ken,
>>
>> Did you try the -o flag ("Override projection check (use current
>> location's projection)
>> Assume that the dataset has same projection as the current location")?
>> https://grass.osgeo.org/grass78/manuals/r.external.html
>> As the error message suggests?
>>
>> If you are sure that projections in fact match, a mismatch in CRS
>> description can be overridden that way, without modifying the input file.
>>
>> Alternatively, you can use r.import to reproject while converting to
>> GRASS format or you create a new location...
>>
>> Cheers
>> Stefan
>>
>>
>> -Original Message-
>> From: grass-user On Behalf Of Ken
>> Mankoff
>> Sent: fredag 23. juli 2021 19:33
>> To: GRASS user list
>> Subject: [GRASS-user] Cannot read NetCDF with 'semi_major_axis' or
>> 'inverse_flattening'
>>
>> Hello list,
>>
>> I have been given an updated NetCDF file - everything is the same as the
>> old version except there are two new attributes in the 'mapping' variable:
>>
>> char mapping ;
>> mapping:geoid = "eigen-6c4" ;
>> mapping:grid_mapping_name = "polar_stereographic" ;
>> mapping:latitude_of_projection_origin = 90. ;
>> mapping:standard_parallel = 70. ;
>> mapping:straight_vertical_longitude_from_pole = -45. ;
>> mapping:semi_major_axis = 6378273. ;
>> mapping:inverse_flattening = 298.27940504282 ;
>> mapping:false_easting = 0. ;
>> mapping:false_northing = 0. ;
>> mapping:_Storage = "contiguous" ;
>>
>>
>> Specifically:
>> mapping:semi_major_axis = 6378273. ;
>> mapping:inverse_flattening = 298.27940504282 ;
>>
>> When I try to load this file with:
>>
>> r.external source=netCDF:file.nc:surface output=surface
>>
>> I get this error:
>>
>> WARNING: Datum not recognised by GRASS and no parameters found
>> ERROR: Projection of dataset does not appear to match current location.
>>
>> Location PROJ_INFO is:
>> name: WGS 84 / NSIDC Sea Ice Polar Stereographic North
>> datum: wgs84
>> ellps: wgs84
>> proj: stere
>> lat_0: 90
>> lat_ts: 70
>> lon_0: -45
>> x_0: 0
>> y_0: 0
>> no_defs: defined
>> init: EPSG:3413
>>
>> Dataset PROJ_INFO is:
>> name: unnamed
>> a: 6378273
>> es: 0.00669388300017
>> proj: stere
>> lat_0: 90
>> lat_ts: 70
>> lon_0: -45
>> x_0: 0
>> y_0: 0
>> no_defs: defined
>>
>> ERROR: datum
>>
>> In case
Re: [GRASS-user] Cannot read NetCDF with 'semi_major_axis' or 'inverse_flattening'
Hi Stefan,
Thanks for the hint. I often use '-o'. In this case after discussion with
the data provider it turns out that it was not EPSG:3413 but instead some
custom projection, even though part of the NetCDF file reported EPSG:3413
in the metadata.
Unfortunately r.import does not work. I get the errors below. Given the
0..3..6..9... it appears to import something, but no raster exists when the
command finishes.
r.import input=NetCDF:file.nc:surface output=surface
WARNING: Datum not recognised by GRASS and no parameters found
WARNING: Projection of dataset does not appear to match current location.
Location PROJ_INFO is:
name: WGS 84 / NSIDC Sea Ice Polar Stereographic North
datum: wgs84
ellps: wgs84
proj: stere
lat_0: 90
lat_ts: 70
lon_0: -45
x_0: 0
y_0: 0
no_defs: defined
init: EPSG:3413
Dataset PROJ_INFO is:
name: unnamed
a: 6378273
es: 0.00669388300017
proj: stere
lat_0: 90
lat_ts: 70
lon_0: -45
x_0: 0
y_0: 0
no_defs: defined
ERROR: datum
WARNING: Datum not recognised by GRASS and no parameters found
WARNING: Datum not recognised by GRASS and no parameters found
Importing raster map ...
0..3..6..9..12..15..18..21..24..27..30..33..36..39..42..45..48..51..54..57..60..63..66..69..72..75..78..81..84..87..90..93..96..99..100
ERROR: Unable to open element file for
ERROR: Unable to get reprojected map extent
-k.
On Sat, Jul 24, 2021 at 2:19 PM Stefan Blumentrath <
stefan.blumentr...@nina.no> wrote:
> Hei Ken,
>
> Did you try the -o flag ("Override projection check (use current
> location's projection)
> Assume that the dataset has same projection as the current location")?
> https://grass.osgeo.org/grass78/manuals/r.external.html
> As the error message suggests?
>
> If you are sure that projections in fact match, a mismatch in CRS
> description can be overridden that way, without modifying the input file.
>
> Alternatively, you can use r.import to reproject while converting to GRASS
> format or you create a new location...
>
> Cheers
> Stefan
>
>
> -Original Message-
> From: grass-user On Behalf Of Ken
> Mankoff
> Sent: fredag 23. juli 2021 19:33
> To: GRASS user list
> Subject: [GRASS-user] Cannot read NetCDF with 'semi_major_axis' or
> 'inverse_flattening'
>
> Hello list,
>
> I have been given an updated NetCDF file - everything is the same as the
> old version except there are two new attributes in the 'mapping' variable:
>
> char mapping ;
> mapping:geoid = "eigen-6c4" ;
> mapping:grid_mapping_name = "polar_stereographic" ;
> mapping:latitude_of_projection_origin = 90. ;
> mapping:standard_parallel = 70. ;
> mapping:straight_vertical_longitude_from_pole = -45. ;
> mapping:semi_major_axis = 6378273. ;
> mapping:inverse_flattening = 298.27940504282 ;
> mapping:false_easting = 0. ;
> mapping:false_northing = 0. ;
> mapping:_Storage = "contiguous" ;
>
>
> Specifically:
> mapping:semi_major_axis = 6378273. ;
> mapping:inverse_flattening = 298.27940504282 ;
>
> When I try to load this file with:
>
> r.external source=netCDF:file.nc:surface output=surface
>
> I get this error:
>
> WARNING: Datum not recognised by GRASS and no parameters found
> ERROR: Projection of dataset does not appear to match current location.
>
> Location PROJ_INFO is:
> name: WGS 84 / NSIDC Sea Ice Polar Stereographic North
> datum: wgs84
> ellps: wgs84
> proj: stere
> lat_0: 90
> lat_ts: 70
> lon_0: -45
> x_0: 0
> y_0: 0
> no_defs: defined
> init: EPSG:3413
>
> Dataset PROJ_INFO is:
> name: unnamed
> a: 6378273
> es: 0.00669388300017
> proj: stere
> lat_0: 90
> lat_ts: 70
> lon_0: -45
> x_0: 0
> y_0: 0
> no_defs: defined
>
> ERROR: datum
>
> In case of no significant differences in the projection definitions, use
> the -o flag to ignore them and use current location definition.
> Consider generating a new location from the input dataset using the
> 'location' parameter.
>
>
> If I remove thos two attributes (semi_major_axis and inverse_flattening)
> with:
>
> ncatted -O -h -a semi_major_axis,mapping,d,l,6378273 -a
> inverse_flattening,mapping,d,f,298.27940504282 file.nc foo.nc
>
> And then try to read in the data, it works.
>
> Is there some reason that those two attributes make it difficult to load
> the data in GRASS? Is there a better work-around than removing them? The
> documented projection of the data has not changed, only some of the
> metadata.
>
> Versions:
> gdal 3.0.4
> grass 7.8.2
>
> Th
[GRASS-user] Cannot read NetCDF with 'semi_major_axis' or 'inverse_flattening'
Hello list, I have been given an updated NetCDF file - everything is the same as the old version except there are two new attributes in the 'mapping' variable: char mapping ; mapping:geoid = "eigen-6c4" ; mapping:grid_mapping_name = "polar_stereographic" ; mapping:latitude_of_projection_origin = 90. ; mapping:standard_parallel = 70. ; mapping:straight_vertical_longitude_from_pole = -45. ; mapping:semi_major_axis = 6378273. ; mapping:inverse_flattening = 298.27940504282 ; mapping:false_easting = 0. ; mapping:false_northing = 0. ; mapping:_Storage = "contiguous" ; Specifically: mapping:semi_major_axis = 6378273. ; mapping:inverse_flattening = 298.27940504282 ; When I try to load this file with: r.external source=netCDF:file.nc:surface output=surface I get this error: WARNING: Datum not recognised by GRASS and no parameters found ERROR: Projection of dataset does not appear to match current location. Location PROJ_INFO is: name: WGS 84 / NSIDC Sea Ice Polar Stereographic North datum: wgs84 ellps: wgs84 proj: stere lat_0: 90 lat_ts: 70 lon_0: -45 x_0: 0 y_0: 0 no_defs: defined init: EPSG:3413 Dataset PROJ_INFO is: name: unnamed a: 6378273 es: 0.00669388300017 proj: stere lat_0: 90 lat_ts: 70 lon_0: -45 x_0: 0 y_0: 0 no_defs: defined ERROR: datum In case of no significant differences in the projection definitions, use the -o flag to ignore them and use current location definition. Consider generating a new location from the input dataset using the 'location' parameter. If I remove thos two attributes (semi_major_axis and inverse_flattening) with: ncatted -O -h -a semi_major_axis,mapping,d,l,6378273 -a inverse_flattening,mapping,d,f,298.27940504282 file.nc foo.nc And then try to read in the data, it works. Is there some reason that those two attributes make it difficult to load the data in GRASS? Is there a better work-around than removing them? The documented projection of the data has not changed, only some of the metadata. Versions: gdal 3.0.4 grass 7.8.2 Thank you, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Extracting data from 4D or higher D (not 3D) NetCDF file
Hello GRASS List, I regularly extract raster slices from a 3D NetCDF using the 'band' option to r.external or r.in.gdal. I now have a NetCDF file with a 4D raster (x,y,layer,t). Is there some way to access all x, all y, one layer, and loop through t? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.watershed identify inland watershed
I'm on mobile and can't fully reply, but isn't there an r.lakes or r.fill module (possibly add on) that does this? Generic answer: r.watershed is the wrong tool, but GRASS has support for this. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Mon, Mar 29, 2021, 11:40 Micha Silver wrote: > Hello: > > You might try `r.param.scale`, or even better `r.geomorphons` modules to > identify geomorphology features, then filter out all pixels identified > as pits. > > > r.watershed is purposely designed to overcome depressions, and find flow > routing thru these spots. So I don't think you can use that module to > identify depressions. > > > On 3/27/21 8:49 PM, ming han wrote: > > Hi Everyone > > > > When I do watershed delineation using r.watershed for great salt > > lake watershed. I found r.watershed always tried to assign an outlet > > for a great salt lake, which does actually not exist because it is an > > inland lake and the great salt lake has no watershed outlet at all. > > > > I noticed that there is a depression option. But is there any > > way that r.watershed can automatically identify depressions while > > defining flow accumulation and stream network? > > > > Thanks > > Ming > > > > ___ > > grass-user mailing list > > grass-user@lists.osgeo.org > > https://lists.osgeo.org/mailman/listinfo/grass-user > > -- > Micha Silver > Ben Gurion Univ. > Sde Boker, Remote Sensing Lab > cell: +972-523-665918 > > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user > ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Data science position in glaciology at GEUS, Copenhagen
Ministry of Climate, Energy and Utilities considers diversity to be
an asset and invites anybody interested, regardless of gender, age, religion or
ethnicity to apply for the position.
* Further information
Please contact the Head of Department Signe Bech Andersen on tel. +45 20 71 34
79, e-mail: s...@geus.dk or senior scientist Ken Mankoff on tel. +45 91 33 38
22 e-mail: k...@geus.dk
International applicants can read more about living and working in Denmark on
GEUS’ homepage:
http://eng.geus.dk/about/jobs-and-education/
* Application
Data science position in Glaciology and Climate (hr-manager.net)
https://candidate.hr-manager.net/ApplicationInit.aspx?cid=5001=140434=5
Applications must include a letter of motivation, a CV, as well as relevant
exam papers, and if relevant, contact details of referees (max. 2).
Your application must be with GEUS no later than 12.00 Noon (Danish time), 29
March 2021. You can send your application by clicking on the "Apply for this
job" ("Søg stillingen") button.
Applications received after the deadline will not be considered.
--
Ken Mankoff (pronouns: any)
Senior Scientist
Geological Survey of Denmark and Greenland (GEUS)
Department of Glaciology and Climate
Copenhagen, Denmark
Mobile (DK): +45 91 33 38 22
Mobile (US): +1 814 826 9682
https://www.geus.dk/om-geus/organisation/afdelinger/glaciologi-og-klima
https://eng.geus.dk/about/organisation/departments/glaciology-and-climate
http://kenmankoff.com
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Polar projection
Hi Markus, On 2020-12-31 at 08:26 -08, Markus Neteler wrote... > On Thu, Dec 31, 2020 at 3:06 AM Ken Mankoff wrote: >> d.grid -g 1:0 color=red >> >> does not show latitude lines on most of the graphic (see attached). > > I managed to generate it with > > Preparation of a 10 degree grid [in EPSG:4326]: > v.mkgrid grid=36,18 map=grid_10deg > > [v.proj from 4326 to polar (3995 or 3413)] > > The resulting map looks as attached (hope I didn't forget to copy a > command here). Indeed, d.grid doesn't looks as expected Your command list was helpful. Yes, things work fine when importing a grid generated in EPSG:4326 to EPSG:3995. Oddly, importing to EPSG:3413 has a different issue (included in the same bug report for now). > d.grid 10 -w .. is incomplete. Worth a bug report, IMO. https://github.com/OSGeo/grass/issues/1224 Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Polar projection
Hi Jón, I work regularly in EPSG:3413, but the same issue applies. "d.grid -g 1:0" doesn't make a grid to the pole. Does it for you? I'd like to submit a bug report as per the suggestion from Markus, but am not sure what the bug is if this command works for you when the pole is included in the view. -k. On 2020-12-31 at 09:23 -08, Jón Eiríksson wrote... > This has worked for me: > > (Thu Dec 31 17:22:11 2020) > g.proj -p > -PROJ_INFO- > name : Stereographic > proj : stere > datum : wgs84 > ellps : wgs84 > lat_0 : 90 > lat_ts : 70 > lon_0 : -45 > k : 1 > x_0: 0 > y_0: 0 > no_defs: defined > towgs84: 0.000,0.000,0.000 > -PROJ_EPSG- > epsg : 3413 > -PROJ_UNITS > unit : meter > units : meters > meters : 1 > (Thu Dec 31 17:22:11 2020) Command finished (0 sec) > > Jon > > > > > On 31 Dec 2020, at 16:26, Markus Neteler wrote: > >> Hi Ken, >> >> On Thu, Dec 31, 2020 at 3:06 AM Ken Mankoff wrote: >>> >>> Dear GRASS List, >>> >>> Can someone suggest what setup to use (EPSG code? proj4 code?) to >>> have the projection centered on the N. pole? I'm looking to create >>> a graphic similar to this image in the Raster Gallery: >>> https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgrass.osgeo.org%2F%2Fimages%2Fgallery%2Fraster%2Fday_on_earth_N.pngdata=04%7C01%7Cjeir%40hi.is%7Cc4ce599d0faf40d3408608d8ada8dff2%7C09fa5f0e211846568529677ed8fdbe78%7C0%7C0%7C637450288118244590%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000sdata=1b0EG%2B%2ByT0LjNGOp55l%2FwYiVKSNfAeUvd%2BZS12dCWZg%3Dreserved=0 >>> >>> If I set up a polar projection based on >>> https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fspatialreference.org%2Fref%2Fsr-org%2F8243%2Fdata=04%7C01%7Cjeir%40hi.is%7Cc4ce599d0faf40d3408608d8ada8dff2%7C09fa5f0e211846568529677ed8fdbe78%7C0%7C0%7C637450288118244590%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000sdata=I7u9KAIJ0geqmtTPZv1vUg3kiAtqtASajvhxzCwFmy4%3Dreserved=0 >>> >>> using: >>> >>> grass -c ./G >>> g.proj -c proj4="+proj=stere +lat_0=90 +lat_ts=45 +lon_0=-170 +k=1 >>> +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" >>> >>> Things mostly work, but >>> >>> d.grid -g 1:0 color=red >>> >>> does not show latitude lines on most of the graphic (see attached). >> >> I managed to generate it with >> >> Preparation of a 10 degree grid: >> >> # EPSG:4326 >> GRASS :~ > g.region -dp >> projection: 3 (Latitude-Longitude) >> zone: 0 >> datum: wgs84 >> ellipsoid: wgs84 >> north: 90N >> south: 90S >> west: 180W >> east: 180E >> nsres: 1 >> ewres: 1 >> rows: 180 >> cols: 360 >> cells: 64800 >> GRASS :~ > v.mkgrid grid=36,18 map=grid_10deg >> # download "Natural Earth I with Shaded Relief, Water, and Drainages" >> from >> https://eur02.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.naturalearthdata.com%2Fdownloads%2F10m-raster-data%2F10m-natural-earth-1%2Fdata=04%7C01%7Cjeir%40hi.is%7Cc4ce599d0faf40d3408608d8ada8dff2%7C09fa5f0e211846568529677ed8fdbe78%7C0%7C0%7C637450288118244590%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000sdata=at%2BPzc3vGV%2FmxVhg%2BxltDODIC52rrl%2Bd39ZTtta7uDA%3Dreserved=0 >> GRASS :~ > r.import in=NE1_HR_LC_SR_W_DR.tif >> output=natural_earth_global_landcover >> GRASS :~ > g.region raster=natural_earth_global_landcover.1 >> GRASS :~ > r.composite r=natural_earth_global_landcover.1 >> g=natural_earth_global_landcover.2 b=natural_earth_global_landcover.3 >> out=natural_earth_global_landcover.rgb >> GRASS :~ > exit >> >> # WGS 84 / Arctic Polar Stereographic >> https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fepsg.io%2F3995data=04%7C01%7Cjeir%40hi.is%7Cc4ce599d0faf40d3408608d8ada8dff2%7C09fa5f0e211846568529677ed8fdbe78%7C0%7C0%7C637450288118244590%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000sdata=JwnkjczUnnrbi%2F%2BnQGepiYo4eB02S5%2BQUWSMK3rRFDs%3Dreserved=0 >> grass78 -c epsg:3995 ~/grassdata/arctic_polar_stereographic >> GRASS :~ > g.proj -w >> PROJCS["WGS
Re: [GRASS-user] GRASS and GEE integration
On 2020-12-17 at 23:45 -08, Bernardo Santos via grass-user wrote... > Hello everyone, I (finally) recently found out about some of the great > capabilities of Google Earth Engine.It is great to be able to get > access to so many satellite images and datasets and process them > without having to download them, build databases, and process them on > local computers. However, I was wondering it would be great to be able > to use all the capabilities from GRASS GIS modules and addons > integrated with GEE.Does anyone know about initiatives to do some kind > of integration in this direction? Thanks!Bernardo There is a Python interface to GEE. There is a Python interface to GRASS. Is that a sufficient bridge? Note that GEE is a Google product, so the chances of it existing in 5 years is relatively slim. 5 years is not long for GRASS. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to find elevation drop for all cells in a basin
I agree this would work. Given that I have 65036 basins, I was trying to avoid a loop. But if it's fast, who cares... I will try this method. Thanks, -k. On 2020-12-16 at 08:29 -08, Thomas Adams wrote... > Hi Ken! > > So, presumably, you know the locations of all the basin outlets > (lat-long, e.g.). Using *v.what.rast*, query the elevation map to get > the elevation at the basin outlet (locations) that the vector points > identify. This associates the basin outlet elevation with its > location. Assuming you have a single elevation map, you'll need to > mask the elevation map to use the method I suggested. You can get the > basin boundary -- if these are not already identified -- to create a > mask using *r.water.outlet*. For each basin, set the MASK and do the > r.mapcalc calculation I proposed, because you know the outlet > elevation and the elevation at all cells within the MASKed basin area. > r.mapcalc will create a new map 'diff' comprised of the elevation > differences between the elevation at the basin outlet and all other > cells in the basin. The elevation difference at the basin outlet will > be ZERO. > > The whole process is easily scriptable. I have done this kind of thing > a fair amount, FWIW... > > Later, after all the diff maps have been created, you could (if you > wanted/needed to) patch them together or do statistics on them, etc. I > hope all this makes sense... > > Best, Tom > > On Tue, Dec 15, 2020 at 3:52 PM Ken Mankoff wrote: > >> >> On 2020-12-15 at 12:28 -08, Thomas Adams wrote... >> > So, if I understand, that seems pretty simple and could be done >> > with r.mapcalc: >> > >> > diff = elev - outlet_elev >> > >> > where elev is a raster map of elevations and outlet_elev is the >> > elevation at the basin outlet. Using a basin mask the calculation >> > could be confined to any basin of interest. But, I am probably >> > missing something, I think… >> >> One of us is missing something, but I'm happy to assume it is me :). >> >> If outlet_elev is only defined at the outlets, diff is only defined >> at the outlets. I want it defined for every inland cell in every >> basin. I have 1,000s of basins. >> >> How do I define outlet_elev everywhere? I can do it now using >> r.stream.distance but prior to Anna suggesting that, I didn't know >> how to do it. Is there some trick like for each basin I set the >> category to the elevation of the outlet, and then use mapcalc @math >> function to create the outlet_elev raster? >> >> -k. >> ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to find elevation drop for all cells in a basin
On 2020-12-15 at 12:28 -08, Thomas Adams wrote... > So, if I understand, that seems pretty simple and could be done with > r.mapcalc: > > diff = elev - outlet_elev > > where elev is a raster map of elevations and outlet_elev is the elevation > at the basin outlet. Using a basin mask the calculation could be confined > to any basin of interest. But, I am probably missing something, I think… One of us is missing something, but I'm happy to assume it is me :). If outlet_elev is only defined at the outlets, diff is only defined at the outlets. I want it defined for every inland cell in every basin. I have 1,000s of basins. How do I define outlet_elev everywhere? I can do it now using r.stream.distance but prior to Anna suggesting that, I didn't know how to do it. Is there some trick like for each basin I set the category to the elevation of the outlet, and then use mapcalc @math function to create the outlet_elev raster? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to find elevation drop for all cells in a basin
Hi Thomas, On 2020-12-15 at 11:34 -08, Thomas Adams wrote... > What exactly do you mean by elevation drop (at each cell)? I can only guess… I meant the change in each cell between their elevation and the elevation of the associated hydrologic outlet cell. I actually need to query a different raster than change in elevation, but it seems like I can pass any cost surface I want as the elevation parameter to r.stream.distance and it works. I think it uses the input direction raster to calculate the outlet for each source, then queries the input and output "elevation" raster (in my case some other cost surface), regardless of what those values represent. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to find elevation drop for all cells in a basin
Hi Anna, On 2020-12-15 at 11:29 -08, Anna Petrášová wrote... > isn't r.stream.distance what you need? Yes that is what I need. I did not realize it worked for the non-stream cells too. Thank you, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] How to find elevation drop for all cells in a basin
Specifically, I have outlets and basins. I can find the (x,y) location of the
outlet for every cell in each basin using the following code. Note that I'm
doing this for all cells in the basin, not just stream cells. The
r.stream.order does provide elev_drop. I'd like this for non-stream cells too,
for each basin.
Assuming that I have a outlets with unique IDs, I can create basins for each
outlet with:
# dir comes from r.stream.extract
r.stream.basins -m direction=dir points=outlets basins=basins
g.copy basins,cat_x
g.copy basins,cat_y
I can then export the x and y location of each outlet to a file:
r.out.xyz input=outlets | awk -F'|' '{print $3, $1}' > ./tmp/outlets_x
r.out.xyz input=outlets | awk -F'|' '{print $3, $2}' > ./tmp/outlets_y
I can then copy my basins (same categories as outlets) and change the category
to x and y.
r.category map=cat_x rules=./tmp/outlets_x separator=space
r.category map=cat_y rules=./tmp/outlets_y separator=space
Note: There is probably a way to combine the above r.out.xyz and r.category
commands so that I don't need to use a temporary disk file.
Finally, convert x and y from category to value
r.mapcalc "outlet_x = @cat_x"
r.mapcalc "outlet_y = @cat_y"
I can then use outlet_x and outlet_y, for example in r.mapcalc to find the
direct distance between each inland cell and the outlet:
r.mapcalc "dist = (outlet_x^2 + outlet_y^2)^(0.5)"
Rather than finding the distance from each inland cell to its outlet, I'd like
to find the elevation change from each inland cell to its outlet. Is this
possible to do efficiently in the Bash interface to GRASS? I can do it with
some loops in Python but prefer to stay in Bash if possible.
Thanks,
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] g.region -a: what's it supposed to do?
Hello again, Please disregard the previous message. I've read the docs again and "even multiples of the resolution" is the key part. I get it now, and why I'm seeing the behavior from the commands below. -k. On 2020-11-27 at 04:13 -08, Ken Mankoff wrote... > My previous understanding is that if I set e,w,n,s AND rows,cols AND > res, and they don't all agree, then res takes precedence if -a is set, > and otherwise res is subject to the other six paramaters. > Alternatively, if just setting e,w,n,s and res, and the e,w res does > not match the n,s res, then setting -a would adjust things. But when > everything agrees, why is -a shifting the bounds? > > grass -c EPSG:3413 G_tmp > > g.region w=-638956 e=856044 s=-3354596 n=-655596 res=1000 -p > eval $(g.region -g) > echo $(( ($n - $s)/1000 )) > echo $(( ($e - $w)/1000 )) > > # why does this shift things? > g.region w=-638956 e=856044 s=-3354596 n=-655596 res=1000 -pa > > > I don't see previous questions about this when searching "g.region resolution > align bounds site:lists.osgeo.org/pipermail/grass-user/" on Google. > > > Thanks, > >-k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] g.region -a: what's it supposed to do?
Hello, I just got caught by a surprising g.region behavior. Can someone clarify what the "-a" flag is supposed to do? From the docs, "Align region to resolution (default = align to bounds, works only for 2D resolution)". I have a defined region on a 1000 m grid. The two are defined together and do "align" as seen by the first set of commands (here "align" means the bounds and the resolution match in both ns and ew directions, and there are no fractional values). In this scenario, why is the application of the "-a" flag shifting the boundaries? What is being aligned to what? My previous understanding is that if I set e,w,n,s AND rows,cols AND res, and they don't all agree, then res takes precedence if -a is set, and otherwise res is subject to the other six paramaters. Alternatively, if just setting e,w,n,s and res, and the e,w res does not match the n,s res, then setting -a would adjust things. But when everything agrees, why is -a shifting the bounds? grass -c EPSG:3413 G_tmp g.region w=-638956 e=856044 s=-3354596 n=-655596 res=1000 -p eval $(g.region -g) echo $(( ($n - $s)/1000 )) echo $(( ($e - $w)/1000 )) # why does this shift things? g.region w=-638956 e=856044 s=-3354596 n=-655596 res=1000 -pa I don't see previous questions about this when searching "g.region resolution align bounds site:lists.osgeo.org/pipermail/grass-user/" on Google. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] [GRASS-dev] Bugs in r.stream.extract
On 2020-11-25 at 04:17 -08, ming han wrote... > And another problem I got is that the flow accumulation I got from > r.accumulate and r.watershed is different when r.accumulate using flow > direction from r.watershed. Again is there anyway r.watershed supports > using flow direction, so we can get the consistent result. There may be reasons for these differences? For example, if SFD v. MFD, or the "-a" flag to r.watershed? > We need this when we need to adjust the flow direction from > r.watershed or r.stream.extract. and then we need to determine new > flow accumulation with an adjusted flow direction dataset. If the > result is inconsistent, not sure what is the solution is. It isn't clear why you are adjusting the flow direction. Is this required? >> But, If I only want to use flow direction to drive streams, which >> function I should use? https://grass.osgeo.org/grass78/manuals/r.water.outlet.html but this only works for 1 outlet. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] [GRASS-dev] Bugs in r.stream.extract
Hi Ming, On 2020-11-23 at 09:05 -08, ming han wrote... > Hope this email finds you well. I got a weird result when using > r.stream.extract, as shown in the following figure. The back grids is > the flow accumulation layer with a flow accumulation threshold larger > than 1000. while the blue line is the stream generated by > r.stream.extract. Why the stream from r.stream.extract did not follow > flow accumulation results? And how to fix this problem? Is there any chance you can share the raster that includes this region so I can examine it? And the exact command you ran? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Reprojecting to x,y grid
On 2020-11-19 at 01:18 -08, Moritz Lennert wrote... > Maybe because of changes linked to proj6 support ? Markus Metz will > know. Yes I found this yesterday. Here is the v.proj recent commits https://github.com/OSGeo/grass/commits/master/vector/v.proj And I'm guessing it was this commit that introduced the behavior change https://github.com/OSGeo/grass/commit/ce56f30a6871356dac0db49c5faa850941e3#diff-fa03e53897bb203999fbd1b78e611058b754555626b19744b8c6f71f95c0955a -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Reprojecting to x,y grid
Hello, On 2020-11-17 at 10:06 -08, Ken Mankoff wrote... > I've successfully worked in rotated pole coordinates (e.g., [1]) but > because I don't know what (x,y) coord contains the pole, I haven't > figured out how to adapt this to the rotated pole situation as in [1]. > Anyway, this isn't a rotated pole... I was incorrect - this is a rotated pole coordinate system. When I run the following code on GRASS 7.4, I get the attached graphic. The key part here is that the rotated pole is at (lon,lat) = (-200,18). grass -c EPSG:4326 Gnorm # Import something into the normal location v.import input=~/data/Zwally_2012/sectors output=Z cat << EOF > ./Gnorm/PERMANENT/PROJ_INFO name: General Oblique Transformation datum: wgs84 towgs84: 0.000,0.000,0.000 proj: ob_tran o_proj: latlon ellps: wgs84 a: 6378137.00 es: 0.0066943800 f: 298.2572235630 lat_0: 0.00 lon_0: 180.00 o_lat_p: 18.0 o_lon_p: -200.0 EOF # rotated_pole:grid_north_pole_latitude = 18. ; # rotated_pole:grid_north_pole_longitude = -200. cat << EOF > ./Gnorm/PERMANENT/PROJ_UNITS unit: degree units: degrees meters: .0174532925 EOF grass -e -c EPSG:4326 Grot grass ./Grot/PERMANENT v.proj location=Gnorm input=Z g.region vector=Z d.mon start=wx0 d.vect Z d.grid 1:0:0 The above code run in GRASS 7.4.0 on a clean Ubuntu 18.04 VM produces the attached graphic. However, when I run the exact same code in GRASS 7.8.4 it does not work. The reproject/rotate does not occur, and Greenland remains at 62-82 N, and 22-72 W. The contents of the "Gnorm" folders is appears identical. The md5sum of DEFAULT_WIND, MYNAME, PROJ_EPSG, PROJ_INFO, and PROJ_UNITS are all the same. I can do the work in GRASS 7.4 in a VM but I'd prefer to do it on the main development machine in the latest GRASS. Does anyone have any idea why this feature stopped working? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Reprojecting to x,y grid
Hi Markus, On 2020-11-17 at 14:40 -08, Markus Neteler wrote... > On Tue, Nov 17, 2020 at 7:06 PM Ken Mankoff wrote: >> I have a data set provided on an (x,y) grid, where each coordinate >> does have a paired (lon,lat). I'd like to work on the native (x,y) >> grid, and bring some other geo-referenced data there. How do I do >> this? > > Here 2.5 wild guesses: > > If it is about "pushing around" a raster map (no squeezing): > - https://grass.osgeo.org/grass78/manuals/r.region.html > > If a rotation is needed: > - https://grass.osgeo.org/grass7/manuals/addons/i.rotate.html > > If it is about a real transform: > - convert the raster matrix to its vector center point representation; then > - https://grass.osgeo.org/grass78/manuals/v.transform.html > - rasterize again Thanks for getting me thinking down this path. This isn't about the raster. I want to work in the raster XY location. I want to reproject a vector there. I know the (lon,lat) of every (x,y). I think maybe I need to solve this outside of GRASS - I can take the (lon,lat) of each point of the geo-referenced vector, and look up the nearest (lon,lat) in the XY location, and set those as the new coordinates for the vector. I'm not yet sure what the difference is between that method and using v.rectify. I would build the control points using the same method: Take the (lon,lat) of each point of the geo-referenced vector, look up the nearest (lon,lat) in the XY location, from there jump to the (x,y) coord in the XY location, and use those as my control points. I then rectify (un-project?) the vector to the XY coords. Is the only difference that the second method would "fill in" points that I don't map between (lon,lat) and (x,y), while the first method requires that I update every point? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Reprojecting to x,y grid
Hi GRASS list, I have a data set provided on an (x,y) grid, where each coordinate does have a paired (lon,lat). I'd like to work on the native (x,y) grid, and bring some other geo-referenced data there. How do I do this? I've successfully worked in rotated pole coordinates (e.g., [1]) but because I don't know what (x,y) coord contains the pole, I haven't figured out how to adapt this to the rotated pole situation as in [1]. Anyway, this isn't a rotated pole... The coordinates for the data look like this: # x, y, lon, lat 1,1,-46.930,59.105 2,1,-46.843,59.127 3,1,-46.757,59.148 4,1,-46.670,59.169 5,1,-46.584,59.191 6,1,-46.497,59.212 7,1,-46.410,59.233 293,441,-10.886,84.128 294,441,-10.415,84.120 295,441,-9.944,84.112 296,441,-9.476,84.104 297,441,-9.008,84.095 298,441,-8.542,84.086 299,441,-8.077,84.076 300,441,-7.614,84.066 Meaning there are 300 columns and 441 rows, each with a (lon,lat). I can easily import this into any projection I want with m.proj, but I'd like to work in this XY coord: grass -c XY ./G g.region s=0 n=441 w=0 e=300 res=1 -pas The problem is I have one other vector file (in EPSG:4326 or EPSG:3413) that I need to use with this grid-coordinate data. How do I reproject or import that? Thanks, -k. [1] https://lists.osgeo.org/pipermail/grass-user/2011-October/062180.html ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Classify basins as "narrow"
On 2020-09-03 at 09:27 -07, Ken Mankoff wrote... > I have skeletons: > > r.to.vect -v input=basins output=basins type=area > v.to.lines input=basins output=bounds Sorry - that's exoskeleton (boundary), not skeleton. The Voronoi skeleton is a powerful tool to simplify the algorithm. Thanks again for this suggestion. > But the voronoi centerline, and then dist from center line to edges, > is a much tighter algorithm than the current one (distance from edges > -> invert -> stream route -> extract stream as centerline -> dist from > center to edges). -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Classify basins as "narrow"
On 2020-09-02 at 23:59 -07, Moritz Lennert wrote... > There is also v.voronoi with the -s flag for extraction of skeletons > from vector polygons. I have skeletons: r.to.vect -v input=basins output=basins type=area v.to.lines input=basins output=bounds But the voronoi centerline, and then dist from center line to edges, is a much tighter algorithm than the current one (distance from edges -> invert -> stream route -> extract stream as centerline -> dist from center to edges). Thank you! -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Classify basins as "narrow"
Hi Stefan, Thank you for these suggestions. I think I need to stick with my more complex workflow for the following reason: Narrowness needs to be computed from the centerline to the edge, rather than from the edge to the center, so that the case of a narrow basin next to a wide basin is properly handled. If you compute from the edge, then the edge of the wide basin would also be considered narrow. Other advantages of the more complex workflow is that the "r.neighbors maximum" suggestion appears to be a binary classifier, while with the more complex method values are actual basin widths, which is better for the downstream fuzzy logic application than a crisp value. Finally, reviewing my code it looks like the "must do this for each basin individually" part at the bottom is incorrect, and the code can be run 1x for all basins. Thank you again for the suggestions, -k. On 2020-09-02 at 23:54 -07, Stefan Blumentrath wrote... > Hei Ken, > > What about a combination of r.grow.distance and r.neighbors? > > 1) Extract boundaries from raster basins > r.neighbors input=basins output=basin_diversity method=diversity > r.reclass input=basin_diversity output=basin_boundary rc=- > "1 = NULL > 2 thru 999 = 1" (if you need the boundaries later on the vector solution > is probably more efficient) > > 2) compute distance from boundaries > r.grow.distance input=basin_boundary output=basin_boundary_distance > > 3) get the maximum distance from the boundaries > r.neighbors method=maximum size="twice what you would consider narrow" > > Depends if it is computationally more efficient... If narrow still > means "a lot" of pixels in width, then r.neighbors might be relatively > slow. But you might increase resolution in that case. > > In order to get the skeleton of the basin you could use r.neighbors as > well: r.neighbors input=basin_boundary_distance method=maximum > output=basin_boundary_distance_max r.mapcalc > expression="skeletons=if(basin_boundary_distance < > basin_boundary_distance_max, null(), 1)" This should approximate the > center lines... Note that you could skip r.neighbors and use the > neighborhood modifier in r.mapcalc... > > Just some ideas for the width issue... > > Cheers, > Stefan ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Classify basins as "narrow"
Hi Markus, On 2020-09-02 at 09:28 -07, Markus Neteler wrote... > On Wed, Sep 2, 2020 at 1:27 PM Ken Mankoff wrote: >> I'd like to detect "narrow" features in GRASS. The attached >> screenshot shows some basins (thick) and streams (thin) and some >> regions (hatched). These regions are spurious because the basin is >> narrow here. I'd like to estimate narrowness with an algorithm. > > Perhaps the compact measure of > > https://grass.osgeo.org/grass78/manuals/v.to.db.html --> compact: > compactness of an area, calculated as compactness = perimeter / (2 * > sqrt(PI * area)) > > could help here? Thank you for sharing this. I did not see your email when I wrote my reply with my solution. I need spatial variance in the result, not a single measure per-basin (that is, I want to be able to say "look out for results in this part of the basin, it is narrow, but don't worry elsewhere in this basin"). However, elsewhere on this project I need to know something about neighboring basins, and was looking into that. Your suggestion to read the v.to.db page led me to the "sides" option there, which I think solves that issue. Thank you! -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Classify basins as "narrow"
Hi List,
On 2020-09-02 at 04:27 -07, Ken Mankoff wrote...
> I'd like to detect "narrow" features in GRASS. The attached screenshot
> shows some basins (thick) and streams (thin) and some regions
> (hatched). These regions are spurious because the basin is narrow
> here. I'd like to estimate narrowness with an algorithm.
>
> I've looked into r.grow.distance r.distance and v.distance but haven't
> been able to imagine a solution yet. Can anyone on this list suggest
> something?
I've solved this, although it takes a lot of steps, and it only computes the
width of the downstream region, not everywhere. The downstream/outflow region
is what I'm interested in (but didn't specify this in my initial post) so this
is OK with me. By this I mean that for a "Y" shaped catchment with three narrow
regions and outflow at the bottom of the stem, the current algorithm will
compute the width of one of the branches and the stem, but not the other branch.
I'm sharing this in case it helps someone else, or someone sees a way to do
this more efficiently (perhaps all basins at once?). Here is the algorithm:
# elevation raster is input
# compute streams, outlets, and basins
r.stream.extract elevation=elevation threshold=11 memory=16384 direction=dir
stream_raster=streams stream_vector=streams
r.mapcalc "outlets = if(dir < 0, 1, null())" # outlets
r.to.vect input=outlets output=outlets type=point
r.stream.basins -m direction=dir points=outlets basins=basins memory=16384 #
basins
# compute distance from edge of each basin
r.to.vect -v input=basins output=basins type=area
v.to.lines input=basins output=bounds
v.to.rast -d input=bounds output=bounds use=val val=1
r.grow.distance input=bounds distance=edge_dist metric=euclidean
# Now the algorithm can only work on one basin at a time.
# Pick a basin with narrow outlet region manually, and mask to that basin.
basin_id=1174
r.mask raster=basins maskcats=${basin_id} --o
# Invert so center line is low, but make outlet lowest.
# Edge dist is 0 at edges, so this cost surface is a hole
# Make the outlet the lowest point so it pours from there.
r.mapcalc "cost_surf = if(isnull(outlets), -edge_dist, -max(edge_dist)-100)"
# Route a stream down the center line
r.stream.extract elevation=cost_surf threshold=11 memory=16384
direction=cost_dir stream_raster=cost_streams stream_vector=cost_streams
# Find the main stream (hack = 1)
r.stream.order -m stream_rast=cost_streams direction=cost_dir
elevation=cost_surf hack=hack memory=16384
r.mapcalc "hack_1 = if(hack == 1, 1, null())"
# Grow center stream to edges
r.grow.distance -m input=hack_1 distance=center_dist
# DONE
Now the value of center_dist at the basin boundary is the distance to the
center line, or 1/2 the width.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.fuzzy.system help with defuzzification
Hi GRASS List,
Following up on the previous message regarding fuzzy logic and GRASS - The
scikit-fuzzy toolbox https://pythonhosted.org/scikit-fuzzy/ is very easy to use
and it looks like it supports a super-set of the functionality provided by the
GRASS addon. The GRASS addon does appear to have some issues, but thanks to
grass-session and the general ease of getting GRASS data into Numpy arrays,
it's easy to work with scikit-fuzzy instead.
Pro's to scikit-fuzzy: maintained, more complete, simpler interface, better
documentation.
Cons: Must leave GRASS, Orders of magnitude slower.
-k.
On 2020-08-31 at 12:12 -07, Ken Mankoff wrote...
> Hello,
>
> https://grass.osgeo.org/grass78/manuals/addons/r.fuzzy.system.html
>
> I'd like to use the =r.fuzzy.system= addons for GRASS, but am having
> trouble understanding the output. I'm concerned there is a bug. Below
> is a minimum example I've created to demonstrate the issue. If anyone
> on this list has experience with these modules or interest in helping
> getting their documentation improved so that they are easier to use, I
> can take the lead but would benefit from some help with this.
>
> I have already improved some of the incorrect docs for these addons
> (see https://github.com/OSGeo/grass-addons/pull/175 ) but am now
> returning to a bug I opened a few months ago here
> https://github.com/OSGeo/grass-addons/issues/176
>
> I'm not 100 % certain there is a bug - I could be using or
> interpreting the results incorrectly. I hope that is the case and this
> is a simple fix. If code changes are required, I will try to implement
> them, but would be grateful for pointers about what/where if someone
> else sees something I have not yet noticed. I am open to including
> collaborators as co-authors on a paper I'm writing if the effort
> needed to get this working and the contribution are significant.
>
> Here (and the GitHub issue above) is the best way I can demonstrate
> the issue:
>
>
> * GRASS Fuzzy Code Test
>
> + Four-cell raster map foo = [1, 2, 3, 4]
> + Foo is "low" when 1, "hi" when 4, and linear change between +
> foo_lo = [1, 0.66, 0.33, 0] + foo_hi = [0, 0.33, 0.66, 1]
> + Output is "sad" when foo is "lo", "happy" when foo is "hi", and
> linear between
> + Expected output:
> + out values = [1, 0.66, 0.33, 0] + out label: [sad=1 happy=0;
> sad=0.66 happy=0.33; sad=0.33 happy=0.66; sad=0 happy=1]
>
> ** Summary
>
> Results from #+BEGIN_SRC bash r.fuzzy.system -m maps=fuzzy.maps
> rules=fuzzy.rules output=out res=10 r.stats -1l out #+END_SRC
>
> Are:
>
> #+RESULTS:
> | 0.63636363 | sad=1.00 | happy=1.00 | 0.54545456 | sad=1.00 |
> | happy=1.00 | 0.36363637 | sad=1.00 | happy=1.00 | 0.27272728 |
> | sad=1.00 | happy=1.00 |
>
> I can see that I've encoded sadness, presumably because that was
> listed first in the fuzzy.maps __OUTPUT__ (that file contents is shown
> below), but both sad=1.0 and happy=1.0 seems incorrect to me.
>
>
> ** Set up domain and helper functions
>
> #+BEGIN_SRC bash :results verbatim grass -c ./G_fuzzy
>
> export GRASS_OVERWRITE=1 export GRASS_VERBOSE=0 g.region w=0 s=0 e=4
> n=1 res=1 -pa -s #+END_SRC
>
> ** Define data
>
> #+BEGIN_SRC bash r.mapcalc "foo = float(col())" r.stats -1 foo
> #+END_SRC
>
> #+RESULTS:
> | 1 | 2 | 3 | 4 |
>
>
> ** Set up rules
>
> #+BEGIN_SRC conf :tangle fuzzy.maps % foo $ lo {right; 1,4; linear; 0;
> 1} $ hi {left; 1,4; linear; 0; 1}
>
> % _OUTPUT_ $ sad {left; 0,1; linear; 0;1} $ happy {right; 0,1; linear;
> 0;1} #+END_SRC
>
> #+BEGIN_SRC conf :tangle fuzzy.rules $ happy { foo = hi } $ sad { foo
> = lo } #+END_SRC
>
> ** Run the fuzzy system with pixel-level debug output
>
> #+BEGIN_SRC bash :results verbatim r.fuzzy.system -m maps=fuzzy.maps
> rules=fuzzy.rules output=out res=10 coors=0.5,0.5 echo " "
> r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10x
> coors=3.5,0.5 #+END_SRC
>
> #+RESULTS: #+begin_example ANTECEDENT happy: 0.000 ANTECEDENT sad:
> 1.000 RESULT (defuzzified): 0.636 UNIVERSE,happy,sad,AGREGATE
> 0.000,0.000,0.000,0.000 0.091,0.000,0.091,0.091
> 0.182,0.000,0.182,0.182 0.273,0.000,0.273,0.273
> 0.364,0.000,0.364,0.364 0.455,0.000,0.455,0.455
> 0.545,0.000,0.545,0.545 0.636,0.000,0.636,0.636
> 0.727,0.000,0.727,0.727 0.818,0.000,0.818,0.818
> 0.909,0.000,0.909,0.909
>
> ANTECEDENT happy: 1.000 ANTECEDENT sad: 0.000 RESULT (defuzzified):
> 0.273 UNIVERSE,happy,sad,AGREGATE 0.000,1.000,0.000,1.000
> 0.091,0.909,0.000,0.909 0.182,0.818,0.000,0.818
> 0.273,0.727,0.000,0.727 0.364,0.636,0.000,0.636
> 0.455,0.545,0.000,0.545 0.545,0.455,0.000,0.455
>
[GRASS-user] Classify basins as "narrow"
Hi GRASS list, I'd like to detect "narrow" features in GRASS. The attached screenshot shows some basins (thick) and streams (thin) and some regions (hatched). These regions are spurious because the basin is narrow here. I'd like to estimate narrowness with an algorithm. I've looked into r.grow.distance r.distance and v.distance but haven't been able to imagine a solution yet. Can anyone on this list suggest something? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.fuzzy.system help with defuzzification
Hello,
https://grass.osgeo.org/grass78/manuals/addons/r.fuzzy.system.html
I'd like to use the =r.fuzzy.system= addons for GRASS, but am having trouble
understanding the output. I'm concerned there is a bug. Below is a minimum
example I've created to demonstrate the issue. If anyone on this list has
experience with these modules or interest in helping getting their
documentation improved so that they are easier to use, I can take the lead but
would benefit from some help with this.
I have already improved some of the incorrect docs for these addons (see
https://github.com/OSGeo/grass-addons/pull/175 ) but am now returning to a bug
I opened a few months ago here https://github.com/OSGeo/grass-addons/issues/176
I'm not 100 % certain there is a bug - I could be using or interpreting the
results incorrectly. I hope that is the case and this is a simple fix. If code
changes are required, I will try to implement them, but would be grateful for
pointers about what/where if someone else sees something I have not yet
noticed. I am open to including collaborators as co-authors on a paper I'm
writing if the effort needed to get this working and the contribution are
significant.
Here (and the GitHub issue above) is the best way I can demonstrate the issue:
* GRASS Fuzzy Code Test
+ Four-cell raster map foo = [1, 2, 3, 4]
+ Foo is "low" when 1, "hi" when 4, and linear change between
+ foo_lo = [1, 0.66, 0.33, 0]
+ foo_hi = [0, 0.33, 0.66, 1]
+ Output is "sad" when foo is "lo", "happy" when foo is "hi", and linear between
+ Expected output:
+ out values = [1, 0.66, 0.33, 0]
+ out label: [sad=1 happy=0; sad=0.66 happy=0.33; sad=0.33 happy=0.66;
sad=0 happy=1]
** Summary
Results from
#+BEGIN_SRC bash
r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10
r.stats -1l out
#+END_SRC
Are:
#+RESULTS:
| 0.63636363 | sad=1.00 | happy=1.00 |
| 0.54545456 | sad=1.00 | happy=1.00 |
| 0.36363637 | sad=1.00 | happy=1.00 |
| 0.27272728 | sad=1.00 | happy=1.00 |
I can see that I've encoded sadness, presumably because that was listed first
in the fuzzy.maps __OUTPUT__ (that file contents is shown below), but both
sad=1.0 and happy=1.0 seems incorrect to me.
** Set up domain and helper functions
#+BEGIN_SRC bash :results verbatim
grass -c ./G_fuzzy
export GRASS_OVERWRITE=1
export GRASS_VERBOSE=0
g.region w=0 s=0 e=4 n=1 res=1 -pa -s
#+END_SRC
** Define data
#+BEGIN_SRC bash
r.mapcalc "foo = float(col())"
r.stats -1 foo
#+END_SRC
#+RESULTS:
| 1 |
| 2 |
| 3 |
| 4 |
** Set up rules
#+BEGIN_SRC conf :tangle fuzzy.maps
% foo
$ lo {right; 1,4; linear; 0; 1}
$ hi {left; 1,4; linear; 0; 1}
% _OUTPUT_
$ sad {left; 0,1; linear; 0;1}
$ happy {right; 0,1; linear; 0;1}
#+END_SRC
#+BEGIN_SRC conf :tangle fuzzy.rules
$ happy { foo = hi }
$ sad { foo = lo }
#+END_SRC
** Run the fuzzy system with pixel-level debug output
#+BEGIN_SRC bash :results verbatim
r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10
coors=0.5,0.5
echo " "
r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10x
coors=3.5,0.5
#+END_SRC
#+RESULTS:
#+begin_example
ANTECEDENT happy: 0.000
ANTECEDENT sad: 1.000
RESULT (defuzzified): 0.636
UNIVERSE,happy,sad,AGREGATE
0.000,0.000,0.000,0.000
0.091,0.000,0.091,0.091
0.182,0.000,0.182,0.182
0.273,0.000,0.273,0.273
0.364,0.000,0.364,0.364
0.455,0.000,0.455,0.455
0.545,0.000,0.545,0.545
0.636,0.000,0.636,0.636
0.727,0.000,0.727,0.727
0.818,0.000,0.818,0.818
0.909,0.000,0.909,0.909
ANTECEDENT happy: 1.000
ANTECEDENT sad: 0.000
RESULT (defuzzified): 0.273
UNIVERSE,happy,sad,AGREGATE
0.000,1.000,0.000,1.000
0.091,0.909,0.000,0.909
0.182,0.818,0.000,0.818
0.273,0.727,0.000,0.727
0.364,0.636,0.000,0.636
0.455,0.545,0.000,0.545
0.545,0.455,0.000,0.455
0.636,0.364,0.000,0.364
0.727,0.273,0.000,0.273
0.818,0.182,0.000,0.182
0.909,0.091,0.000,0.091
#+end_example
** Run the fuzzy system to get raster output
#+BEGIN_SRC bash
r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10
#+END_SRC
#+RESULTS:
Test the intermediate result from =-m= flag.
#+BEGIN_SRC bash
r.stats -1 out_sad
#+END_SRC
#+RESULTS:
| 1 |
| 0.6669 |
| 0.3334 |
| 0 |
#+BEGIN_SRC bash
r.stats -1 out_happy
#+END_SRC
#+RESULTS:
| 0 |
| 0.3334 |
| 0.6669 |
| 1 |
Now inspect the final result
#+BEGIN_SRC bash :results verbatim
r.stats -1l out
#+END_SRC
#+RESULTS:
: 0.63636363 sad=1.00 happy=1.00
: 0.54545456 sad=1.00 happy=1.00
: 0.36363637 sad=1.00 happy=1.00
: 0.27272728 sad=1.00 happy=1.00
+ The intermediate results, =out_sad= and =out_happy= make sense, but the final
result does not.
+ This suggests there is an issue with the defuzzification step. Let's try
other methods.
#+BEGIN_SRC bash :results verbatim
r.fuzzy.system -m maps=fuzzy.maps rules=fuzzy.rules output=out res=10
defuz=bisector imp=minimum
r.stats -1l out
#+END_SRC
#+RESULTS:
:
: 0.63636363 sad=1.00 happy=1.00
: 0.54545456
Re: [GRASS-user] r.surf.area
Hi Valter, On Fri, May 29, 2020 at 3:43 AM Valter Albino wrote: > Suppose you compute the "r.surf.area"(1) to some raster file and find that > the 3d area is smaller than your 2d area? > What could be the problem? > Note: > The exercise was done with a 5 m cell size raster file, with 37 km2 > watershed with a slope of 0,0341 m/m, within 7.8 version of GRASS GIS 7.8 > console in a 3.12.3 QGIS [with this method the difference is higher (12 ha, > comparing with ArcGIS)] and outside QGIS, in GRASS GIS GUI > Thanks in advance > Can you provide a small example, perhaps just a few grid cells, that reproduces the problem? Also, it isn't clear if you see the same problem in ArcGIS, QGIS, and GRASS, or just GRASS. How are you calculating the 2D area? Do you note in r.surf.area under DESCRIPTION that it says " Therefore, area of a flat surface will be reported as (rows + cols -1) * (area of cell) less than area of flat region due to a half row and half column missing around the perimeter." and does this explain what you're seeing? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Global Hydrological modeling with r.watershed, r.stream.extract, r.stream.extract
Hi Guiseppe, I've successfully run with 4.5 billion cells. How many cells do you have? I notice you do not have the "-m" flag to tell it to use disk swap in place of all memory. Maybe that would help? -k. On 2020-05-05 at 11:55 -07, Giuseppe Amatulli wrote... > Dear GRASS Team > > I am running a global analysis where I need to use "tiles" as computational > units in which I use the following three commands: > r.watershed -b elevation=elv depression=dep accumulation=flow > drainage=dir_rw flow=pixel_area memory=10 --o --verbose > r.stream.extract elevation=elv accumulation=flow depression=dep > threshold=0.05 direction=dir_rs stream_raster=stream memory=10 --o > --verbose > r.stream.basins -l stream_rast=stream direction=dir_rs basins=lbasin > memory=10 --o --verbose > > The basins that were not completely within a tile (resulting in > broken-basins) have been removed (see below the three tiles in Figs. 1,2,3 > including only entire basins), > and now I'm in the phase of merging all the tiles having only complete > basins. > > When I merge the tiles (Fig 4), some basin borders do not match perfectly, > and some areas have NoData (see Fig 5,6) or have the Basin ID of the below > basin (Fig 7). > I noticed that these phenomena appear only when I merge tiles that have > very large broken-basins that can not be included in the tile due to RAM > limitations. > > My thought is that r.stream.basins needs the entire dimension of two > adjacent basins to be able to detect the border without gap and without a > potential random selection. > Is there any part of the r.stream.basins code that I can potentially check > and eventually hack to avoid this problem? > > For the rest, all the RAM limitation and other problem have been solved > soon we will have a global stream network and basin delineation performed > 100% in GRASS!!! > > Thank you > Best Regards > Giuseppe > > > > Fig 1. Left Tile > [image: image.png] > > > Fig 2. Center tile > [image: image.png] > > > Fig 3. Right Tile > [image: image.png] > > > Fig 4. Merge all the tiles > [image: image.png] > > > Fig 5. Gap -> small white area > [image: image.png] > > > Fig 6. Gap -> small white area > [image: image.png] > > Fig 7. boarder Basins inconsistency among tiles > > [image: image.png] ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] garray.read() returns 0-filled array
Hello,
Replying to my earlier issues from one month ago:
On 2020-03-15 at 08:36 -07, Ken Mankoff wrote...
> foo = garray.array().read("foo", null=np.nan)
I have found I can read the GRASS array into Python with this syntax:
foo = garray.array(mapname="foo", null=np.nan)
I'm still waiting to upgrade GRASS to test Stefan's suggestion (later/elsewhere
in this thread).
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] garray.read() returns 0-filled array
Hi Stefan, On 2020-03-16 at 07:01 -07, Stefan Blumentrath wrote... > Hei Ken, > > With this PR: > https://github.com/zarch/grass-session/pull/14 > applied, also pygrass modules should work as usual with grass-session. Thank you for submitting that. It'll take me a while to get up and running with the git version as opposed to the pip version, but I will test it. > GRASS GIS 7.8 should use Python 3, so I am surprised to see that your > $GRASS_PYTHON version is /usr/bin/python2.7 ... That is possibly also > a source for your trouble with garray... I'm on an LTS OS - Ubuntu 18.04 and the default GRASS that comes with that: 7.4. I'm looking forward to 7.8 but haven't upgraded yet. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] garray.read() returns 0-filled array
On 2020-03-15 at 15:11 -07, Stefan Blumentrath wrote... > I have no experience with garray, but raster2numpy function in pygrass > worked quite well for me: > https://grass.osgeo.org/grass78/manuals/libpython/pygrass.raster.html#pygrass.raster.raster2numpy Thank you for this suggestion. This works but raises new issues. 1) Fully external-to-GRASS interactions works better for me. This solution requires running Python from within GRASS, which adds an extra layer, which makes it harder to get to the same place interactively in Python and debug. 2) Worse, it seems I need to run the $GRASS_PYTHON version which is /usr/bin/python2.7, so I no longer have access to my scientific Python stack. The 1st issue can be worked around. The 2nd not so much. I hope I'm mis-using the library and can do this from within my Anaconda Python. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] garray.read() returns 0-filled array
Hello,
I'm trying to use GRASS via Python and import some rasters into numpy. I'm
using garray.read(raster) and the result is always 0. When I inspect the raster
in GRASS (via bash) it is not zero. I can re-create this with the following
MWE. Can someone suggest what I might be doing wrong or if this is a bug?
Thank you,
-k.
import numpy as np
import shutil
import os
from grass_session import Session
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
import grass.script.array as garray
if os.path.exists("G_test"): # reset everything
shutil.rmtree("G_test", ignore_errors=False, onerror=None)
W = Session()
W.open(gisdb='.', location='G_test', create_opts="EPSG:3413")
g.region(w=0, s=0, n=100, e=100)
r.mapcalc("foo = 42.24")
foo = garray.array().read("foo", null=np.nan)
W.close()
print(np.max(foo))
# 0.0 is
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Can I turn off compression
Hi,
Thanks for the replies. Yes I found LZ4 to be the best. In answer to Vaclav
suggesting general interest, here are some code and results so you can test
this on your own system which will have different disk/CPU/RAM speeds. I am
only reporting time, not disk usage here.
I'm not totally sure how to find the size of a GRASS raster from within GRASS.
I could go to the OS to find it. A better report would include a scatter plot
of time v. size, but I only cared about time.
Summary: Use LZ4 compression within GRASS. For any maps that already exist, it
isn't clear if LZ4 or uncompressed is faster. They seem about equal on my
system using this as a test:
r.mapcalc "compressed = 42"
r.mapcalc "uncompressed = 42"
r.compress -u map=uncompressed
/usr/bin/time -f %E parallel -j 1 -N 0 r.mapcalc '"foo = compressed * 2"' :::
$(seq 10)
/usr/bin/time -f %E parallel -j 1 -N 0 r.mapcalc '"foo = uncompressed * 2"' :::
$(seq 10)
-k.
grass -c EPSG:3413 ./G
g.region w=0 e=1 s=0 n=1 res=1 -pa
export GRASS_OVERWRITE=1
for C in RLE LZ4 BZIP2; do
export GRASS_COMPRESSOR=${C}
echo -n "${C} ":
/usr/bin/time -f %E parallel -j 1 -N 0 r.mapcalc '"foo = 42"' ::: $(seq 10)
done
export GRASS_COMPRESSOR=ZLIB
for ZL in $(seq -1 10); do
export GRASS_ZLIB_LEVEL=${ZL}
echo -n "ZLIB ${ZL} ":
/usr/bin/time -f %E parallel -j 1 -N 0 r.mapcalc '"foo = 42"' ::: $(seq 10)
done
RLE :0:17.39
LZ4 :0:14.51
BZIP2 :0:24.05
ZLIB -1 :0:22.49
ZLIB 0 :0:19.14
ZLIB 1 :0:20.31
ZLIB 2 :0:20.05
ZLIB 3 :0:19.87
ZLIB 4 :0:22.72
ZLIB 5 :0:22.49
ZLIB 6 :0:22.25
ZLIB 7 :0:22.32
ZLIB 8 :0:22.26
ZLIB 9 :0:22.67
ZLIB 10 :0:19.77
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Can I turn off compression
Hm. Answer: No I just did: export GRASS_COMPRESSOR=none And then get this error message: WARNING: No compression is not supported for GRASS raster maps, using default ZLIB Seems like a shame to not have this as an option. I'll try speeding it up by tweaking the compression method and amount. -k. On 2020-02-29 at 08:37 -08, Ken Mankoff wrote... > Hi, > > I'd like to turn off compression for all maps. Is this possible with > g.gisenv or environment variables? The g.gisenv manual only shows how > to choose among compression methods, not disable it. > > This does not seem to work: > > g.gisenv set="GRASS_COMPRESSOR=none" # optimize for speed > > I can decompress one map with: > > r.compress -u map=foo > > But I'm trying to skip the time spent compressing and decompressing > altogether. I could use a small speed boost, and this seems like an > easy way to speed things up - I can afford the disk space for all the > temporary intermediate maps. > > Thanks, > > -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Can I turn off compression
Hi, I'd like to turn off compression for all maps. Is this possible with g.gisenv or environment variables? The g.gisenv manual only shows how to choose among compression methods, not disable it. This does not seem to work: g.gisenv set="GRASS_COMPRESSOR=none" # optimize for speed I can decompress one map with: r.compress -u map=foo But I'm trying to skip the time spent compressing and decompressing altogether. I could use a small speed boost, and this seems like an easy way to speed things up - I can afford the disk space for all the temporary intermediate maps. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] compute path with a maximum slope
r.cost generates paths. Set high cost to traverse cells with high slopes? -k. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Tue, Dec 10, 2019, 10:41 li...@lazygranch.com wrote: > I have a DEM in grass using a proj that is in UTM. Viewshed, slope, etc > work fine. What I want to do is have grass compute a path between two > points such that the slope doesn't exceed a maximum value. I don't care > about the path length. > > I suspect there is some guide on the internet for this problem, but > just can't compose the right question for a search engine. A link would > be appreciated. If not, some basic guidance would be helpful. > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Import CORDEX data in GRASS
Hi Johannes, On 2019-10-23 at 10:49 -07, Johannes Radinger wrote... > So maybe another way (if rotated pole is not possible with GRASS/gdal) Rotated pole works fine in GRASS. Here are my notes from when I did this a few years ago. -k. https://lists.osgeo.org/pipermail/grass-user/2011-October/062180.html #+BEGIN_SRC bash :results verbatim # Normal and rotated GRASS locations trash Gnorm Grot grass74 -c EPSG:4326 Gnorm # Import something into the normal location v.import input=~/data/Zwally_2012/sectors output=Z cat << EOF > ./Gnorm/PERMANENT/PROJ_INFO name: General Oblique Transformation datum: wgs84 towgs84: 0.000,0.000,0.000 proj: ob_tran o_proj: latlon ellps: wgs84 a: 6378137.00 es: 0.0066943800 f: 298.2572235630 lat_0: 0.00 lon_0: 180.00 o_lat_p: 18.0 o_lon_p: -200.0 EOF # rotated_pole:grid_north_pole_latitude = 18. ; # rotated_pole:grid_north_pole_longitude = -200. cat << EOF > ./Gnorm/PERMANENT/PROJ_UNITS unit: degree units: degrees meters: .0174532925 EOF grass74 -e -c EPSG:4326 Grot grass74 ./Grot/PERMANENT v.proj location=Gnorm input=Z g.region vector=Z d.mon start=wx0 d.erase d.vect Z d.grid 1:0:0 #+END_SRC ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing .jpg & georectifying
Hi Rich, A maybe-easier workflow is to load QGIS with an XYZ Google Satellite base layer, then save to image *with a World File*. That should allow easy import into other tools. https://geogeek.xyz/how-to-add-google-maps-layers-in-qgis-3.html On 2019-10-18 at 12:52 -07, Rich Shepard wrote... > I have a .jpg screenshot of a Google-Earth Pro view and see that r.in.gdal > supports that file format. However, I have no idea how to structure the > command line for the task. > > On the GE view I recorded the geographic coordinates of four > readily-identifiable points and can put them in a v.in.ascii point file for > import > > My web searches for the process of saving a GE view as an .kml file were > without success. That would be another option for importing the GE image and > using it as the background for additional layers outlining specific areas. > > Pointers to appropriate modules is appreciated. > > Regards, > > Rich > > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Maps in a defined region
On 2019-10-03 at 15:25 +02, Ondřej Pešek wrote...
> but r.out.gdal would export "something" from all the maps in the loop,
> wouldn't it? I thought that if I would loop through all my rasters and
> pass them to r.out.gdal, it would create an output for every of them,
> just the ones outside the computational region will be full of no data
> values. I wanted to avoid this export of empty rasters.
In the loop you cant test each raster for some non-null values:
for r in $(g.list...); do
eval $(r.univar -g ...)
if [[ ${sum} != 0 ]]; then...
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Basin DEM advice needed
Hi Rich, On 2019-10-03 at 00:22 +02, Rich Shepard wrote... > Attached are two maps using 1m LiDAR data. The annotated map, > basin-elevations.png, was drawn by individually applying d.rast to each of > the 70 maps covering the basin. Sharp breaks can be seen where the data > cross quads or flights didn't match up smoothly. > > The un-annotated map, nehalem-dem-patched.png, displays the results of > running r.patch on all 70 maps. Topographically it's quite different from > the individual maps; almost flat when the north, east, and south edges > should have elevations similar to the other map. The two maps are showing the same data. d.rast scales each raster to (min,max) of the colorbar, which is why the tiles appear misaligned. The data is unchanged. > I think I should apply r.resamp.stats to aggregate the 1m resolution > to 5m. I'd like your thoughts on this. I don't see how anyone can answer that without knowing the specifics of your project and goals. > I assume that I should resample each individual map, then re-run > r.patch on the coarser maps because r.slope.aspect and r.info need a > single map as input. I'd patch than resample to avoid possible edge effects, if your resample isn't perfectly aligned with each sub-raster. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.patch does not find list of raster maps
On 2019-10-01 at 18:55 -04, Rich Shepard wrote... > # get the list of DEM raster files > MAPS = `g.list type=raster sep=newline pat="45123*"` MAPS=$(g.list type=raster sep=comma pat="45123*") I think $() is preferred over `` and you need a comma separated list for r.patch. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Resample DEM to larger cells
On 2019-10-01 at 17:48 -04, Rich Shepard wrote... > There are 70 individual raster maps within the mask of the basin's > watershed boundary. Am I correct that speed-issues aside, I need to > run r.buildvrt so I have a single map as input to modules such as > r.info (to identify the highest elevation), r.slope, and others? While > I understand that a mask limits analyses to its area that is based on > a single map extending beyond the mask and not a collection of 70 > 'tiles'. Am I correct? r.univar can take multiple multiple maps as input and from that you could get max. r.slope and others would work best on a VRT or r.patch from the individual rasters. > Are there alternatives to a virtual raster for situations like this? r.patch? Or for loop over all the individual rasters w/ r.info to find the max? That won't work for r.slope because of edge effects. A VRT is probably easiest / best. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Importing DEM data with different projections
On 2019-09-17 at 14:41 -07, Rich Shepard wrote... > My question [...] is how do I handle these situations? Do I just > override the dataset's projection and force it to use the existing > location's? Overriding only makes sense to me if the data *is* in the same projection but the metadata is incorrect. > I created a separate location for each map then reprojected each to my > statewide location. I'm sure there's a more efficient way to handle > these discrepancies. r.import would save you creating the 2nd map projection and a few lines of code. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] qgis and grass prosseses error.
Hi Håvard,
The first thing I would do is to try to remove non-ASCII characters and spaces
from your path. Paths like this:
r.external input="C:\Gis\øving 8\klippet transformert.tif"
And the å in your name may be causing problems.
-k.
On 2019-09-13 at 07:33 -07, havard solheim wrote...
> Hello im student named Håvard taking a gis class in Norway. I keep
> getting the same error when trying to use grass in qgis and my
> teatcher told me to ask here since he had no Clue.
>
> Under is the log from trying to do a V.rast.stats in qgis. The program
> Version are also here im not sure what to post either in this mail
> since i am really New to this.
>
>
>
> QGIS version: 3.8.1-Zanzibar
>
> QGIS code revision: dcd95cc648
>
> Qt version: 5.11.2
>
> GDAL version: 2.4.1
>
> GEOS version: 3.7.2-CAPI-1.11.0 b55d2125
>
> PROJ version: Rel. 5.2.0, September 15th, 2018
>
> Processing algorithm…
>
> Algorithm 'v.rast.stats' starting…
>
> Input parameters:
>
> { 'GRASS_MIN_AREA_PARAMETER' : 0.0001, 'GRASS_OUTPUT_TYPE_PARAMETER' : 0,
> 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None,
> 'GRASS_SNAP_TOLERANCE_PARAMETER' : -1, 'GRASS_VECTOR_DSCO' : '',
> 'GRASS_VECTOR_EXPORT_NOCAT' : False, 'GRASS_VECTOR_LCO' : '', 'column_prefix'
> : 'hogst', 'map' : 'C:/Gis/øving 8/hedmark_kommuner.gpkg|layername=output',
> 'method' : [0,1,2,3,4,5,6,7,8,9,10,11,12], 'output' : 'TEMPORARY_OUTPUT',
> 'percentile' : 90, 'raster' : 'C:/Gis/øving 8/klippet transformert.tif' }
>
>
>
> g.proj -c proj4="+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
> +units=m +no_defs"
>
> v.in.ogr min_area=0.0001 snap=-1.0 input="C:\Gis\øving
> 8\hedmark_kommuner.gpkg" layer="output" output="vector_5d7ba884cdc462"
> --overwrite -o
>
> r.external input="C:\Gis\øving 8\klippet transformert.tif" band=1
> output="rast_5d7ba884cde683" --overwrite -o
>
> g.region n=6957956.983371275 s=6633986.983371275 e=717076.6256832004
> w=529936.6256832004 res=30.0
>
> v.rast.stats -c map=vector_5d7ba884cdc462 raster=rast_5d7ba884cde683
> column_prefix="hogst"
> method="number,minimum,maximum,range,average,stddev,variance,coeff_var,sum,first_quartile,median,third_quartile,percentile"
> percentile=90 --overwrite
>
> v.out.ogr type="auto" input="vector_5d7ba884cdc462"
> output="C:/Users/Håvard/AppData/Local/Temp/processing_7868b90be3504757845921771ddee21e/aee2bf1504c74baaa7b62cf0308f635e/output.gpkg"
> format="GPKG" --overwrite
>
> Cleaning up temporary files...
>
> Starting GRASS GIS...
>
> WARNING: Concurrent mapset locking is not supported on Windows
>
> Executing
>
> ...
>
> C:\Users\H†vard\Documents>chcp 1252 1>NUL
>
> C:\Users\Håvard\Documents>g.proj -c proj4="+proj=utm +zone=32 +ellps=GRS80
> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
>
> Default region was updated to the new projection, but if you have multiple
> mapsets `g.region -d` should be run in each to update the region from the
> default
>
> Projection information updated
>
> C:\Users\Håvard\Documents>v.in.ogr min_area=0.0001 snap=-1.0
> input="C:\Gis\øving 8\hedmark_kommuner.gpkg" layer="output"
> output="vector_5d7ba884cdc462" --overwrite -o
>
> Over-riding projection check
>
> Check if OGR layer contains polygons...
>
> 0..4..9..13..18..22..27..31..36..40..45..50..54..59..63..68..72..77..81..86..90..95..100
>
> Creating attribute table for layer ...
>
> WARNING: Unable to open database
>
> by driver
>
> ERROR: Unable to open database
>
> by driver
>
> C:\Users\Håvard\Documents>r.external input="C:\Gis\øving 8\klippet
> transformert.tif" band=1 output="rast_5d7ba884cde683" --overwrite -o
>
> Over-riding projection check
>
> Reading band 1 of 1...
>
> Link to raster map created.
>
> C:\Users\Håvard\Documents>g.region n=6957956.983371275 s=6633986.983371275
> e=717076.6256832004 w=529936.6256832004 res=30.0
>
> C:\Users\Håvard\Documents>v.rast.stats -c map=vector_5d7ba884cdc462
> raster=rast_5d7ba884cde683 column_prefix="hogst"
> method="number,minimum,maximum,range,average,stddev,variance,coeff_var,sum,first_quartile,median,third_quartile,percentile"
> percentile=90 --overwrite
>
> WARNING: Coor file of vector map is larger
> than it should be (18 bytes excess)
>
> WARNING: Coor file of vector map is larger
> than it should be (18 bytes excess)
>
> Topology not available for vector map .
> Registering primitives...
>
> WARNING: Unable to open vector map on level
> 2. Try to rebuild vector topology with v.build.
>
> ERROR: Unable to open vector map
>
> ERROR: An error occurred while converting vector to raster
>
> WARNING: No data base element files found
>
> Execution of
>
> finished.
>
> Cleaning up temporary files...
>
> Press any key to continue . . .
>
> Cleaning up temporary files...
>
> Starting GRASS GIS...
>
> WARNING: Concurrent mapset locking is not supported on Windows
>
> Executing
>
> ...
>
> C:\Users\H†vard\Documents>chcp 1252 1>NUL
>
> C:\Users\Håvard\Documents>v.out.ogr type="auto" input="vector_5d7ba884cdc462"
[GRASS-user] resample STRDS to finer temporal resolution
Hello, I'd like to take an STRDS with several overlapping time span rasters and extract a time stamp (or a shorter time-span) that is some operation of the overlapping rasters. For example: % t.rast.list input=SEC | head name|mapset|start_time|end_time SEC_19920101_19970101|CCI_SEC|1992-01-01 00:00:00|1997-01-01 00:00:00 SEC_19930101_19980101|CCI_SEC|1993-01-01 00:00:00|1998-01-01 00:00:00 SEC_19940101_19990101|CCI_SEC|1994-01-01 00:00:00|1999-01-01 00:00:00 SEC_19950101_2101|CCI_SEC|1995-01-01 00:00:00|2000-01-01 00:00:00 SEC_19960101_20010101|CCI_SEC|1996-01-01 00:00:00|2001-01-01 00:00:00 And I'd like to get a raster that represents some subset, for example a time stamp of 1994-07-01 or time span covering that day, or a time span covering one month. The operation could be simple - say the average of all rasters that include or contain the timestamp, or I'd be happy if I could have more complex operations, like a weighted running mean. I've looked at t.rast.aggregate - this appears to only work for sampling to larger time spans. For example, I can aggregate by decade: % t.rast.aggregate -n input=SEC output=SEC_y basename=base granularity="10 years" method=average % t.rast.list input=SEC_y name|mapset|start_time|end_time base_1992|CCI_SEC|1992-01-01 00:00:00|2002-01-01 00:00:00 base_2002|CCI_SEC|2002-01-01 00:00:00|2012-01-01 00:00:00 base_2012|CCI_SEC|2012-01-01 00:00:00|2022-01-01 00:00:00 But not by day: % t.rast.aggregate -n input=SEC output=SEC_d basename=base granularity="10 days" method=average % t.rast.list input=SEC_d ERROR: Space time raster dataset not found If I want this (copied from above): % t.rast.list input=SEC | head name|mapset|start_time|end_time SEC_19920101_19970101|CCI_SEC|1992-01-01 00:00:00|1997-01-01 00:00:00 SEC_19930101_19980101|CCI_SEC|1993-01-01 00:00:00|1998-01-01 00:00:00 SEC_19940101_19990101|CCI_SEC|1994-01-01 00:00:00|1999-01-01 00:00:00 SEC_19950101_2101|CCI_SEC|1995-01-01 00:00:00|2000-01-01 00:00:00 SEC_19960101_20010101|CCI_SEC|1996-01-01 00:00:00|2001-01-01 00:00:00 resampled monthly for a decade, do I need to do it all manually? For example, loop over 120 months, and for each one 1) t.rast.extract or t.rast.mapcalc to get STRDS rasters that contain that month, 2) r.mapcalc to generate a new raster where I use t.rast.list as the input to the expression, and then 3) add the new raster to a new monthly STRDS? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Creating location using EPSG code
On 2019-08-27 at 16:42 +02, Rich Shepard wrote... > While I know that regions are subsets within the bounds of the CRS > defined for a location I don't recall importing a vector or raster map > and not being able to zoom to its extents before explicitly setting a > region. Perhaps in the earlier work the region was set w/o import based on how you created the location? If you created the location using "-c SomeGeoTiff.tif", then the region is set to that file. If you create it with "-c EPSG:", it doesn't know extent, resolution, etc. and the defaults, as you saw, are all 0 or 1. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Copy projection to another location
On 2019-08-26 at 22:38 +02, Rich Shepard wrote... > can I change the CRS for an existing location or only when it's being > created? You can, but it is rare that you should do it, and a fairly low-level operation. It is required in some edge cases, like when setting up a rotated pole coordinate system. What is your use-case or intended goal/behavior/setup? Perhaps there are easier and more appropriate methods. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Creating location using EPSG code
Hi Rich, On 2019-08-26 at 23:10 +02, Rich Shepard wrote... > The region in the project location is incorrect. No. Setting the region and location are different things. You set the location with your "-c EPSG:2838". > With only one small data set imported I decided to start over by > deleting that location in /data/grassdata/. You haven't shown the import command, but you can set the region from imported data with: g.region (raster=... | vector=) Or explicitly: g.region e=... w=... n=... s=... res=... If setting with vector=... you'll need to set the "res" explicitly. > What have I done incorrectly? Not set your region, just your location. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] wms google
On 2019-08-05 at 12:14 -04, Sebastián Dietrich wrote... > is it possible to load xyz tiles of google images as in qgis? Yes. See "No plugin required" answer, not accepted answer, here: https://gis.stackexchange.com/questions/20191/adding-basemaps-from-google-or-bing-in-qgis > Did somebody get it work in grass? I haven't tried yet. I did just try exporting/saving a Google XYZ layer from QGIS to my Desktop as a projected GeoTIFF and it did not seem to work. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Fwd: Re: remove small islands from vector areas
On 2019-06-21 at 21:10 +02, Robert Nuske wrote...
> If there is no way to only treat the islands/holes (which i don't know
> about), it will change the outer boundary as well. Fjord like
> structures will vanish. Neighboring areas might even merge if close
> enough. And since one can not control the buffer style of sharp edges
> (like in postgis, geos), everything will become "rounder".
Can you work in raster space or do you have to remain with resolution-agnostic
vectors? I've solved this issue with maintaining coastlines in raster space by
finding and saving the largest (or several largest) clumps (e.g. the ocean
surrounding my domain), filling all holes which may involve filling the sea or
growing the coastline, and then re-applying the saved clumps, thereby removing
any artificially rounded coast.
r.clump input=raster output=clumps
clump_ID_largest=$(r.stats -c clumps sort=desc | head -n1 | cut -d" " -f1)
r.mapcalc "mask_ocean = if(clumps == ${clump_ID_largest}, null(), 1)"
r.mfilter -z input=raster output=raster_overfilled filter=filter.txt
r.mapcalc "raster_filled = raster_overfilled * mask_ocean" --o
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] too many categories: buffer overflow
On 2019-06-20 at 09:37 +02, Moritz Lennert wrote... > Just a rapid, wild guess here, but would it be feasible to just > vectorize the three separately and then create the link afterwards > using v.distance ? This may work. An alternate method is if I do this through the Python GRASS interface, I can more easily create unique IDs that are not based on the simple "Cell #" algorithm I'm currently using. That is harder to do in bash. Or just recompile GRASS to use LONGs for primary key rather than INT. For now I'm sticking with a lower resolution 90 m raster. Problem solved. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] too many categories: buffer overflow
Hi Micha and Markus,
On 2019-06-18 at 10:07 -04, Micha Silver wrote...
> Do you really want a vector polygon map with > 2 billion features?
No, and there are not that many.
% r.info -r basins
min=-2147474681
max=2147429730
But I don't have categories from 1 to 2147429730. The values are sparse. I
describe my workflow and why I've created these sparse values in more detail
below.
Even though << 2 billion, there should be many basins. This is all of Greenland
at 30 m resolution, which is 4.5 billion features.
Taking a step back, I'm trying to generate unique basin values that match the
stream and outlet CAT values. Here is my workflow which doesn't appear to have
any problems when run at 90x90 m resolution (400 million cells) but fails at
30x30 m resolution (10x as many, or 4.5 billion cells).
1) Find streams:
r.stream.extract elevation=head threshold=${THRESH} memory=16384 direction=dir
stream_raster=streams stream_vector=streams
2) Find outlets. Where streams have outlets, use the same CAT value so the two
can be linked in further analysis. But many outlets don't have streams. These
need to have unique categories for the next step when we find basins. This is
where my error is. I set the unique value to the cell #, which is > 2 billion
when using a 30x30 m domain.
r.mapcalc "outlets_all = if(dir < 0, 1, null())"
r.mapcalc "outlets_streams_1 = if((dir < 0) && (not(isnull(streams))), streams,
outlets_all)"
### BUG INTRODUCED HERE, setting (eventual) cat to cell number:
r.mapcalc "outlets_streams = if(outlets_streams_1 != 1, outlets_streams_1,
max(outlets_streams_1)+1+col()+(max(col())*(row()-1)))"
# convert outlets to a vector.
r.out.xyz input=outlets_streams | \
v.in.ascii input=- output=outlets_streams separator=pipe \
columns="x int, y int, cat int" x=1 y=2 cat=3
Q: How can I create the outlets_streams vector for all locations where dir < 0
(all outlets), that maintains the same value as the streams raster where that
raster is defined, but unique values at all other locations where streams is
not defined, but dir < 0?
3) Find basins
r.stream.basins -m direction=dir points=outlets_streams basins=basins_all
memory=16384 --verbose
4) Absorb small basins
r.clump -d input=basins_all output=basins_nosmall minsize=124
r.mode base=basins_nosmall cover=basins_all output=basins
### BUG APPEARS HERE
r.to.vect -v input=basins output=basins type=area
# drop outlets for absorbed basins.
r.mapcalc "outlets = if(outlets_streams == basins, basins, null())"
r.to.vect -v input=outlets output=outlets type=point
NOTE: I use r.mode instead of r.area because I need to maintain the category
value, so that eventual vectors can have linked primary keys. r.area re-assigns
categories.
Any advice how to generate streams, outlets, and basins all with linked primary
key would be much appreciated.
Thanks,
-k.
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] too many categories: buffer overflow
Hello, I think I'm experiencing a buffer overflow. This is a hard one to search for with GRASS GIS because the word "buffer" and "overflow" appear throughout as in r.buffer and overflowing weirs, etc. but I'm referring to the C-code error type of buffer oveflow: r.to.vect -v input=basins output=basins type=area DBMI-SQLite driver error: Error in sqlite3_step(): UNIQUE constraint failed: basins.cat ERROR: Unable to insert into table: insert into basins values ( -2137121269, '(Category -2137121269)') Does anyone have any suggestion how I can convert a raster with MANY categories to a vector? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] create basins of similar area
Maybe make basins with threshold an order of magnitude smaller. Then merge small areas to get your 4-5 km2 areas. Still not totally sure how this makes sense hydrologically. -k. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Thu, Apr 25, 2019, 16:05 Markus Metz wrote: > > > On Thu, Apr 25, 2019 at 5:34 PM Moritz Lennert < > mlenn...@club.worldonline.be> wrote: > > > > > > > > Am 25. April 2019 17:29:24 MESZ schrieb Martin Landa < > landa.mar...@gmail.com>: > > >Hi, > > > > > >čt 25. 4. 2019 v 17:19 odesílatel Francois Chartier > > > napsal: > > >> Use r.wathershed or r.basins. > > > > > >well, I know these tools. I am just not sure how to control optimal > > >basin area (eg. 4-5km2). I just found `threshold` parameter. > > > > > > > > > I don't really see how a homogeneous size could actually make sense. The > threshold allows you to define a minimum size, but depending on the actual > topography actual sizes of basins will always differ. > > I agree. A solution could be manual editing: merge basins by adhering > hydrological rules (merge a downstream basin with an upstream basin granted > that other upstream basins do not become isolated). But hydrologically this > does not make sense, compare with stream order numbers. > > Markus M > > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
On 2019-03-17 at 23:28 -0700, Moritz Lennert wrote... > Note that you are still using the -c flag, and so at each export all > the geometries are exported, even those that have no category value in > the respective layer. And that was the cause of the issue. The output file no longer appear in triplicate in QGIS or file size. They still compress well. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
On 2019-03-17 at 11:13 -0700, Markus Metz wrote... > [...] GeoPackage file sizes are indeed confusing. In my experience, > compressing GeoPackage files can usually reduce file size > substantially. Yes, that helps a lot. 739 MB to 150 MB. I'll do this. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
On 2019-03-17 at 10:31 -0700, Saber Razmjooei wrote... > You might need to run VACUUM from gdal tools: > > ogrinfo -sql "VACUUM" pal.gpkg That did not change the size. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
On 2019-03-17 at 10:03 -0700, Ken Mankoff wrote... > After learning from Markus & Moritz I see that the GPKG files are as > they should be - vector features are not in triplicate. I think there > are some strange QGIS import issues where it imports everything 3x but > this is not a GRASS issue. Actually, I am still confused about something with the file sizes. The data/metadata all appear correct. But if I export a three-layer vector w/ three tables 1x, the GPGK file is 12 MB and contains 1 table. It is missing two tables. If I then export each of the databases associated with each of the layers, they are 0.2 MB. If export the vector 3x, once per layer, with update, as advised earlier in this thread, the resulting GPKG is ~34 MB, or ~3x as large as I would expect it to be. v.out.ogr -c input=pal output=pal_all.gpkg --o # 12 MB, but only layer1 table db.out.ogr input=pal output=l1.csv layer=1 # 7.5k db.out.ogr input=pal output=l2.csv layer=2 # 124K db.out.ogr input=pal output=l3.csv layer=3 # 85k v.out.ogr -c input=pal output_layer=tbl2_points output=pal.gpkg layer=1 --o v.out.ogr -c input=pal output_layer=tbl2_areas output=pal.gpkg layer=2 -u --o v.out.ogr -c input=pal output_layer=tbl2_lines output=pal.gpkg layer=3 -u --o # pal.gpg is 34 MB, with 3 tables. Why does adding two additional (small) tables increase the size of the resulting geopackage by 100s of times the size of the tables? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
Hello, After learning from Markus & Moritz I see that the GPKG files are as they should be - vector features are not in triplicate. I think there are some strange QGIS import issues where it imports everything 3x but this is not a GRASS issue. Thanks, -k. On 2019-03-16 at 09:37 -0700, Moritz Lennert wrote... > Le Sat, 16 Mar 2019 09:20:17 -0700, Ken Mankoff a > écrit : > >> Hi Moritz [and Markus] >> >> On 2019-03-16 at 08:24 -0700, Moritz Lennert >> wrote... >> > Try running v.out.ogr three times, changing layer each time and >> > using the -u flag for updating the existing gpkg file: >> > >> > v.out.ogr -c input=pal output_layer=pal1 output=pal.gpkg layer=1 >> > v.out.ogr -c input=pal output_layer=pal2 output=pal.gpkg layer=2 -u >> > --o v.out.ogr -c input=pal output_layer=pal3 output=pal.gpkg >> > layer=3 -u --o >> >> That seems to work. When I then read in the GeoPackage and diff, >> things look better. >> >> v.in.ogr input=pal.gpkg output=pal2 --o diff <(v.db.connect -g pal) >> <(v.db.connect -g pal2) # no major diff, but table name is not >> preserved. >> >> But the vector data is also included 3x now. This may be due to >> v.patch? I'm adding Markus because he suggested the method to patch >> different layers. After I patch the three layers: >> >> rm -fR ./nc_spm_08_grass7/gpkg grass ./nc_spm_08_grass7/PERMANENT >> g.mapset -c gpgk g.copy vector=poi_names_wake@PERMANENT,points g.copy >> vector=lakes@PERMANENT,areas g.copy vector=railroads@PERMANENT,lines >> v.category in=areas out=areas_2 op=chlayer layer=1,2 --o v.category >> in=lines out=lines_3 op=chlayer layer=1,3 --o v.patch >> in=points,areas_2,lines_3 output=pal --o >> >> I would expect that v.info reports different information if I use the >> "layer" option. But all three layers report IDENTICAL info. This diff >> command is silent: >> >> diff3 -A <(v.info pal layer=1) <(v.info pal layer=2) <(v.info pal >> layer=3) > > By default v.info provides geometry information which is indeed the > same, as the map contains all these geometries. IIUC, the layer option > for v.info is meant for the v.info -c flag to decide which table's > columns to display. > >> >> Also, the final export commands suggested by Moritz: Layers are >> v.out.ogr -c input=pal output_layer=tbl_pts output=pal.gpkg layer=1 >> --o v.out.ogr -c input=pal output_layer=tbl_areas output=pal.gpkg >> layer=2 -u --o v.out.ogr -c input=pal output_layer=tbl_lines >> output=pal.gpkg layer=3 -u --o >> >> seem to export everything 3x. This is based on 1) the message >> "Exporting 15279 areas (may take some time)..." repeated each time 2) >> the file size increasing by 12 MB each time, and 3) reading in, >> displaying, and clicking brings up everything 3x in the "Query >> Results" window. >> >> Should I just keep my vectors and databases each as a separate vector >> and not use layers? >> > > That is my fault: I shouldn't have used the -c flag. Try without it. > You will see different feature counts in each layer when you run > > ogrinfo pal.gpkg -so pal1 > > Moritz ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] GeoPackage and multiple tables
Hi Moritz [and Markus] On 2019-03-16 at 08:24 -0700, Moritz Lennert wrote... > Try running v.out.ogr three times, changing layer each time and using > the -u flag for updating the existing gpkg file: > > v.out.ogr -c input=pal output_layer=pal1 output=pal.gpkg layer=1 > v.out.ogr -c input=pal output_layer=pal2 output=pal.gpkg layer=2 -u --o > v.out.ogr -c input=pal output_layer=pal3 output=pal.gpkg layer=3 -u --o That seems to work. When I then read in the GeoPackage and diff, things look better. v.in.ogr input=pal.gpkg output=pal2 --o diff <(v.db.connect -g pal) <(v.db.connect -g pal2) # no major diff, but table name is not preserved. But the vector data is also included 3x now. This may be due to v.patch? I'm adding Markus because he suggested the method to patch different layers. After I patch the three layers: rm -fR ./nc_spm_08_grass7/gpkg grass ./nc_spm_08_grass7/PERMANENT g.mapset -c gpgk g.copy vector=poi_names_wake@PERMANENT,points g.copy vector=lakes@PERMANENT,areas g.copy vector=railroads@PERMANENT,lines v.category in=areas out=areas_2 op=chlayer layer=1,2 --o v.category in=lines out=lines_3 op=chlayer layer=1,3 --o v.patch in=points,areas_2,lines_3 output=pal --o I would expect that v.info reports different information if I use the "layer" option. But all three layers report IDENTICAL info. This diff command is silent: diff3 -A <(v.info pal layer=1) <(v.info pal layer=2) <(v.info pal layer=3) Also, the final export commands suggested by Moritz: v.out.ogr -c input=pal output_layer=tbl_pts output=pal.gpkg layer=1 --o v.out.ogr -c input=pal output_layer=tbl_areas output=pal.gpkg layer=2 -u --o v.out.ogr -c input=pal output_layer=tbl_lines output=pal.gpkg layer=3 -u --o seem to export everything 3x. This is based on 1) the message "Exporting 15279 areas (may take some time)..." repeated each time 2) the file size increasing by 12 MB each time, and 3) reading in, displaying, and clicking brings up everything 3x in the "Query Results" window. Should I just keep my vectors and databases each as a separate vector and not use layers? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] GeoPackage and multiple tables
Hello, I'm trying to export a vector to GeoPackage. The vector has 3 layers and each layer has an associated table. When I view the data (pre-export) in GRASS it looks good, but when I export and open in QGIS I only see the layer1 table. Am I doing something wrong or is this a bug/limit in the GeoPackage exporter? Here is an MWE using the nc_spm_08_grass7 data set. I'm using grass 7.4.0. I note v.out.ogr does warn > # WARNING: The combination of types is not supported by all formats. But I think this is always printed and that GeoPackage can support a combination of types. -k. grass ./nc_spm_08_grass7/PERMANENT g.mapset -c gpgk g.copy vector=poi_names_wake@PERMANENT,points g.copy vector=lakes@PERMANENT,areas g.copy vector=railroads@PERMANENT,lines v.category in=areas out=areas_2 op=chlayer layer=1,2 v.category in=lines out=lines_3 op=chlayer layer=1,3 v.patch in=points,areas_2,lines_3 output=pal v.db.addtable map=pal table="tbl_points" layer=1 v.db.addtable map=pal table="tbl_areas" layer=2 v.db.addtable map=pal table="tbl_lines" layer=3 v.db.connect -p map=pal # Vector map is connected by: # layer <1/tbl_points> table in database through driver with key # layer <2/tbl_areas> table in database through driver with key # layer <3/tbl_lines> table in database through driver with key # GRASS 7.4.0 (nc_spm_08_grass7):~/data/grass > v.db.addcolumn map=pal layer=1 column="pt INT" v.db.addcolumn map=pal layer=2 column="ar INT" v.db.addcolumn map=pal layer=3 column="ln INT" v.db.update map=pal layer=1 column=pt value=42 v.db.update map=pal layer=2 column=ar value=43 v.db.update map=pal layer=3 column=ln value=44 v.out.ogr -c input=pal output=pal.gpkg --o # WARNING: The combination of types is not supported by all formats. -- Ken Mankoff Geological Survey of Denmark and Greenland (GEUS) Department of Glaciology and Climate ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Seeking vector layer tutorial
Hi Markus, Thank you for your reply and help. I'm starting to gain a better understanding of the power and complexity of the GRASS vector implementation. -k. On 2019-03-10 at 21:09 +0100, Markus Metz wrote... > On Fri, Mar 8, 2019 at 9:22 PM Ken Mankoff wrote: >> >> Hello, >> >> I'm trying to gain proficiency w/ GRASS vectors & layers. I've read > https://grass.osgeo.org/grass76/manuals/vectorintro.html and > https://grasswiki.osgeo.org/wiki/Vector_Database_Management a few > times, but am still not clear on how to create layers. >> >> I have three vector layers: BASINS, STREAMS, and OUTLETS. All basins >> have > outlets. All streams have outlets. Not all basins have streams (some > small basins have no streams but still have an outlet). >> >> I can merge these three with v.patch, but they're all in one layer. I > thought it might make sense for each to exist on its own layer. Or > perhaps my misconception is that they cannot - it is the attribute > tables that exist on separate layers? >> >> If someone can provide an example how to combine these there vectors, >> or > link to a tutorial, I would be grateful. > > In this example, there are three vectors, each with one layer which is > probably for all layer number 1. > > If you want to patch these three vectors together such that each input > vector with one layer ends up in one individual layer in the output > vector, you need to do some preparation, e.g.: > > Keep vector layer 1 for BASINS. > > Change vector layer 1 for STREAMS to layer 2 with v.category > in=STREAMS out=STREAMS_l2 op=chlayer layer=1,2 > > Change vector layer 1 for OUTLETS to layer 3 with v.category > in=OUTLETS out=OUTLETS_l3 op=chlayer layer=1,3 > > patch them together with v.patch in=BASINS,STREAMS_l2,OUTLETS_l3 > out=... > > In order to copy attributes, you need to copy the attribute tables of > the input vectors BASINS, STREAMS, and OUTLETS. > > The attribute tables linked to the original input vectors are printed > with v.db.connect -p or v.db.connect -g. > > Now copy the tables with db.copy and link them to corresponding layer > of the patched output using v.db.connect. > > HTH, > > Markus M ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] merging small areas but keep raster categories
Hi Stefan, I just solved this with the following and was replying to my own message when yours came in. r.clump -d input=basins output=basins_nosmall minsize=123 r.mode base=basins_nosmall cover=basins output=basins_merged On 2019-03-10 at 12:06 +0100, Stefan Blumentrath wrote... > Did you try the rmarea method for r.reclass.area? That converts to > vector, runs v.clean with rmarea and then converts back to raster (if > I understood correctly... The issue with rmarea is I don't want to remove small areas, I want to merge them with larger ones. If there is an island of only small ones, they are all dropped because there is nothing large to attach them to. The r.clump code I use above correctly merges them into one, and then the r.mode line re-assigns the most common (largest area) original category value. The latter (keeping the main category value) was what I was having trouble with when I wrote my original message. Thank you for thinking about the issue and replying, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] merging small areas but keep raster categories
I'd like to merge small areas from a raster with their large neighbors. I can do this with `r.clump` with the `minsize` argument. But r.clump reclassifies everything, which is causing problems downstream with my workflow. I've tried working around this with `r.reclass.area ... mode=greater`, but that doesn't have the same functionality. Can anyone suggest a workflow or tool that merges clumps but does not re-categorize them? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Seeking vector layer tutorial
Hello, I'm trying to gain proficiency w/ GRASS vectors & layers. I've read https://grass.osgeo.org/grass76/manuals/vectorintro.html and https://grasswiki.osgeo.org/wiki/Vector_Database_Management a few times, but am still not clear on how to create layers. I have three vector layers: BASINS, STREAMS, and OUTLETS. All basins have outlets. All streams have outlets. Not all basins have streams (some small basins have no streams but still have an outlet). I can merge these three with v.patch, but they're all in one layer. I thought it might make sense for each to exist on its own layer. Or perhaps my misconception is that they cannot - it is the attribute tables that exist on separate layers? If someone can provide an example how to combine these there vectors, or link to a tutorial, I would be grateful. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Adding map to layer manager (wxpython GUI) from grass console?
On 2019-02-07 at 13:58 +0100, Joe wrote... > An other question: what way would you suggest to focus on for learning > purpose? Monitors and terminal console, or a GUI based approach? I'll offer opposing view. Terminal. CLI is more powerful than GUI. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] get the shifted minimum from a raster series
Hi Frank, On 2019-02-07 at 08:27 +0100, Frank David wrote... > I try to find how to get the shifted minimum value of a series of > raster map. I want to shift by one, or more, the "minimum" value found > on each cell to get the "almost minimum" of the series. Is there an > available function/method to do that ? I don't think I fully understand what you want. Are you working with strds? You mention "series". If I wanted to get the second-minimum of a single raster, I would use r.stats with the "-1" flag, sort, and search for the 2nd (or 2nd-to-last). -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
Hi Nikos, On 2019-02-03 at 18:45 +0100, Nikos Alexandris wrote... > there is also > https://grass.osgeo.org/grass77/manuals/r.blend.html. Thanks - looks useful. I ended up with: r.mapcalc and "r#..." to split 5 different rasters into RGB channels. r.patch'ing 3x (5x r, 5x g, 5x b) then rgb displaying the patched rasters in ps.map. r.blend looks easier. But I'm not sure how to keep things spread evenly and looking good for >2 rasters. Is there any deep technical reason ps.map does not accept >1 raster? I don't want to start GUI programming, but from the command line it seems like this could be a tractable first major contribution to the GRASS code, depending on the reasons it doesn't exist. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
On 2019-02-02 at 11:51 +0100, Ken Mankoff wrote... > One more ps.map question. What is the best/easiest method to patch > together multiple rasters with different colors? From https://grasswiki.osgeo.org/wiki/Ps.map_scripts """If you want to show two or more overlapping raster maps you need to combine them with the r.patch module or r.mapcalc's '#' color operator. (see also the r.his and r.composite modules)""" -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
Hello, On 2019-01-29 at 11:48 +0100, st_kie...@web.de wrote... > Give it a try. ps.map is not that complicated. The manual page even > offers an out of the box template for instant results (just scroll to > the end). And if you get familiar with that tool you can produce very > convincing graphics for publications, with graphically additions of > your liking. One more ps.map question. What is the best/easiest method to patch together multiple rasters with different colors? I have two raster basemaps plus a velocity raster map. One basemap is the output of r.shade (a colored shaded relief map). One basemap is the output of r.relief (a grayscale shaded map) The velocity raster has its own color scale. I'd like all three in ps.map, but understand there is a one-raster limitation (any good reason for this?). How do I patch them but maintain distinct colors? Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
On 2019-01-29 at 18:42 +0100, Veronica Andreo wrote... > I just tested with the cartographic composer GUI (g.gui.psmap) and > with command line ps.map, and as long as the other mapset is > accessible, no problems with fully qualified map names, i.e., > map@mapset. No need to copy maps, you can use map@mapset in the > instruction file OK - maybe there was something else that was causing the problem for me. Thank you for testing and letting me know this does work. -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
On 2019-01-29 at 11:48 +0100, st_kie...@web.de wrote... > Give it a try. ps.map is not that complicated. The manual page even > offers an out of the box template for instant results (just scroll to > the end). And if you get familiar with that tool you can produce very > convincing graphics for publications, with graphically additions of > your liking. You're right - it is easy. One issue - I can't reference map@mapset. I have to g.copy map@mapset,map and then reference just "map". Am I missing something or is this a real limitation? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to generate a vector graphics file from GRASS CLI?
I will look into ps.map. I was hoping for something simpler. Please excuse brevity. Sent from tiny pocket computer with non-haptic feedback keyboard. On Tue, Jan 29, 2019, 09:28 Hi Ken, > have you considered using ps.map? > > https://grass.osgeo.org/grass76/manuals/ps.map.html > > regards > > Stefan > > > Ken Mankoff hat am 29. Januar 2019 um 08:43 > geschrieben: > > > > > > Hi, > > > > I'm trying to generate a vector graphics file (PDF, but EPS or PS is > fine) from the GRASS CLI. I'd like to have both raster and vector output. > When I follow the simple examples from PSDRIVER or CAIRO the result is > always the same - only the last graphic appears. It is as though there is a > "d.erase" between every command. > > > > From https://grass.osgeo.org/grass76/manuals/psdriver.html > > > > export GRASS_RENDER_IMMEDIATE=ps > > export GRASS_RENDER_TRUECOLOR=TRUE > > g.region raster=elevation > > d.rast elevation > > d.vect roadsmajor color=red > > > > Only shows roads, no elevation. > > > > I can get what I want with the PNG driver, but those graphics are not > publication quality. > > > > Thanks, > > > > -k. > > ___ > > grass-user mailing list > > grass-user@lists.osgeo.org > > https://lists.osgeo.org/mailman/listinfo/grass-user > ___ > grass-user mailing list > grass-user@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] How to generate a vector graphics file from GRASS CLI?
Hi, I'm trying to generate a vector graphics file (PDF, but EPS or PS is fine) from the GRASS CLI. I'd like to have both raster and vector output. When I follow the simple examples from PSDRIVER or CAIRO the result is always the same - only the last graphic appears. It is as though there is a "d.erase" between every command. From https://grass.osgeo.org/grass76/manuals/psdriver.html export GRASS_RENDER_IMMEDIATE=ps export GRASS_RENDER_TRUECOLOR=TRUE g.region raster=elevation d.rast elevation d.vect roadsmajor color=red Only shows roads, no elevation. I can get what I want with the PNG driver, but those graphics are not publication quality. Thanks, -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Overlapping pixels among DEM tiles to compute the LS factor for RUSLE
Hi Nikos, On 2019-01-25 at 07:18 -0800, Nikos Alexandris wrote... > A billion-pixel scaled DEM is the main input to compute the slope length > and steepness (LS) factor for RUSLE (`r.watershed`), only. > > Tiling this DEM in tiles of 5K^2 pixels (`r.tile`), appears to be a > reasonable approach to parallelise this process. Do you need to parallelise it? I just ran r.watershed on a 4.5 billion-pixel DEM w/o a problem on my laptop and I think it took ~6 hours, but I'm not sure. It might have been closer to 12. I have 32 GB of ram and gave half to the process, but told it to use disk instead of memory via the -m flag: r.watershed -s -m elevation=phi threshold=100 drainage=dir memory=16384 --v --o -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] New Google Scholar profile: "GRASS Development Team"
On 2018-12-09 at 10:24 -0800, Markus Neteler wrote: > On Sun, Dec 9, 2018 at 7:20 PM Ken Mankoff wrote: >> On 2018-12-09 at 02:24 -0800, Markus Neteler wrote: >> > for fun, to easily track GRASS GIS related publications I have >> > created a new Google Scholar profile: GRASS Development Team (a >> > virtual team!) >> > https://scholar.google.com/citations?user=gJ0ZB0cJ >> >> Is this papers about GRASS? > > yes > >> Or papers that use GRASS? > > also yes > >> And papers by anyone, > > yes > >> or papers by the GRASS dev team? > > not limited to GRASS dev team since that quite changed in the past 35+ > years. > > So: "openly GRASS GIS related" was the idea. Should that be modified? I'm not involved enough that I think my suggestions should carry weight, but I'll make an observation. The current list is all-inclusive and any of those paper scan be found by searching for "GRASS GIS" on Google Scholar. All this adds is a citation count & h-index. Which is a fine point and may be the point. A perhaps useful curated subset is papers *about* GRASS (or papers where GRASS is the focus, add-ons, etc). That curated list is probably important enough it is worth being hosted as a page on the GRASS wiki, and not just on Google Scholar. Anyway, the current focus includes me: non-dev, just user, used it in a paper, so feel free to add DOI http://dx.doi.org/10.5194/tc-11-303-2017 -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] New Google Scholar profile: "GRASS Development Team"
On 2018-12-09 at 02:24 -0800, Markus Neteler wrote: > for fun, to easily track GRASS GIS related publications I have created > a new Google Scholar profile: GRASS Development Team (a virtual team!) > https://scholar.google.com/citations?user=gJ0ZB0cJ Is this papers about GRASS? Or papers that use GRASS? And papers by anyone, or papers by the GRASS dev team? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] vector segment lengths
On 2018-12-09 at 03:49 -0800, Markus Neteler wrote: > I have made an attempt to add that (r73781): > https://grass.osgeo.org/grass77/manuals/v.split.html#notes-on-segment-length-information > > Pls fix as needed, then I'll backport it. How do you prefer small documentation change suggestions. I have grass_trunk via svn. Post an email here? Patch here? Patch to grass-dev? Issue on Trac? -k. ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
