Re: [GRASS-user] Re: Is it possible to export a Raster to KML
Roger is correct. GE will just drape your raster data over it's terrain mesh regardless of what vertical parameters your data has. It is essentially 2D with a z value and is interpreted as such. That is why EPSG:4326 works, because GE doesn't care about any other information. No conversions/shifts occur. I export rasters from GRASS, generate hillshades, and overlay in GE using gdal2tiles daily. I seldom see any alignment issues, and when I do it is usually my data. Of course there are alignment issues in GE, but they will vary place to place regardless of topography. - Jamie On Fri, Oct 28, 2011 at 2:06 PM, Roger André ran...@gmail.com wrote: I find this thread interesting. In my experience overlaying both raster (orthophoto) and vector data in Google Earth with KML, EPSG:4326 has worked with no problems whatsoever. So far as I know, GE does not apply a vertical datum, it simply drapes the 2D data over its internal terrain model. Would it be possible for you to share some data where you have seen that this is a problem? Thanks, Roger On Fri, Oct 28, 2011 at 1:05 PM, Hamish hamis...@yahoo.com wrote: Helmut wrote: ... in my experience overlaying kml (vector and raster) with EPSG:4326 in google earth works relatively well in flat areas, in mountain areas there are sometime distortions therefore. Markus: I can confirm this from my (limited) experience to export data to KML. To correct for the geoid heights, you could use the previously indicated map and calculate the local ellipsoid-geoid difference and vertically shift your vector data with v.transform. Then export to KML. For raster, use r.mapcalc to shift. n.b., IIUC raster KML support in google earth is essentially a hack building on/abusing their raster-icon support. they may have improved things, but this is what I was lead to believe at the time of writing r.out.kml. probably it is useful to consider if the ellipsoid-geoid difference is important vs. the overall vertical differences in the data, if the result is just for viewing maybe it is not worth the trouble. Hamish ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] wiki link to addon r.bilateral broken
Hello all, I'd like to install the addon r.bilateral from the add-ons page ( http://les-ejk.cz/files/programs/grass/r.bilateral.tgz) but the link is broken. Does anyone have the source tgz they could share with me? Thanks, Jamie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Processing Atlas Lidar Data - Problem(s)
Hi Tim, When dealing with lidar data, I always import using r.in.xyz first. You can easily vary the resolution to generate a quick preview and it will give you a sense of the data source and what information it contains. If you then find you need to do more advanced processing, you can move over to the lidar tools. Cheers, Jamie On Thu, Mar 10, 2011 at 10:23 AM, TimNorwey timmy_wey...@live.at wrote: Hi all! I´m trying to generate a DEM/DSM out of LIDAR data that I have downloaded from the website http://atlas.lsu.edu/lidar/ . The data can be downloaded as .csv-file. Within this file no attribute names are given. For test purpose I manipulated the .csv, so that only 10,000 points were left and I added a row at the top x,y,z as attribute names. After doing that I imported the .csv into QGIS using the Plugin Add Delimited Text Layer, stored the points as .shp-file and imported them into GRASS GIS database using v.in.ogr (as 3D points). Then I tried to follow the workflow described on the website http://grass.osgeo.org/wiki/LIDAR#LIDAR_Tools (I started at DEM/DSM separation the more complex way). The first step should be the use of v.extract to separate first and last pulse laser points. Because of the reason that I have no information about first/last pulse laser points, I skipped this step and tried to follow the other steps. But I think this was for nothing, because v.lidar.edgedetection needs the last pulse laser points and v.lidar.growing needs the first pulse laser points. So I don´t write more about this... I really don´t know how to generate a DEM/DSM based on this data and it would be nice if anyone of you has an idea. Is there a possibility to do generate a DEM/DSM without using the steps that are described within the WIKI page above? What do you think about the import of the .csv into GRASS GIS via QGIS? Is there an easier way to import it? If you need any further information to understand my problem, just let me know. Tim -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Processing-Atlas-Lidar-Data-Problem-s-tp6158817p6158817.html Sent from the Grass - Users mailing list archive at Nabble.com. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] How to merge rasters using feathering/blending ?
I wrote a script to do something like this at one point, though I can't find it now. It was intended for terrain, but could be used per band on imagery. It used r.buffer to create a distance gradient and then r.recode to convert this gradient to a percentage. I then used r.mapcalc to merge the 2 rasters based on the percentage. The script automated the buffer distances and the creation of the recode.txt file, but went something like this: ## make a null raster r.mapcalc temp_null = if(isnull(source),1,null) ## create a distance gradient from the null raster. The output will have int values 1-7 where 1 equals the input null raster r.buffer in=temp_null out=temp_buf distances=10,20,30,40,50,60 units=meters ## fill nulls in the distance gradient with 0 for recoding r.mapcalc temp_buf_fill=if(isnull(temp_buf),0,temp_buf) ## recode to a percentage r.recode in=temp_buf_fill out=temp_recode rules=recode.txt where recode.txt would have the range from temp_buf_fill expressed as percentages: 1:1:100:100 2:2:86:86 3:3:71:71 4:4:57:57 5:5:43:43 6:6:29:29 7:7:14:14 0:0:0:0 ## merge the output r.mapcalc merged=if(isnull(source),patch,(patch*(temp_recode/100.0))+(source*(1.0-(temp_recode/100.0 Sorry if that doesn't work, but hopefully you get the idea. - Jamie On Thu, Feb 24, 2011 at 10:53 AM, kaipi mapcoll...@gmx.net wrote: I need to patch rasters containing NULL areas with alternative rasters. I can do this with r.patch but due to different quality/resolution of the rasters this produces prominent seamlines. Is there a way to feather/blend between merged rasters using Grass GIS? Thanks, Kaipi -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/How-to-merge-rasters-using-feathering-blending-tp6061489p6061489.html Sent from the Grass - Users mailing list archive at Nabble.com. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Opensource and pretty maps design anybody?
I generally open GRASS data in QGIS, create a basic map and export to svg (or pdf), and then do the finishing touches in Inkscape. This workflow has worked very well for me. - Jamie On Tue, Aug 17, 2010 at 8:56 AM, kapo coulibaly kmcou...@gmail.com wrote: I've been using opensource GIS for a while (GRASS a lot). But they are all very limited when it comes to creating nice maps. I always have to resolve to using ArcGIS to put final touches on maps. Anybody has a suggestion for a free software to use for map design? ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] dems from coordinate lists
You want the points to represent cell centers so you need to expand your region 1/2 pixel in all four directions. Also look at using r.in.xyz, it works directly on this type of data. Jamie On Thu, May 13, 2010 at 7:26 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I've obtained DEMs in text files with columns X, Y and Z at 25m spacing. The first three entries in the text file are: - X,Y,Z 99550,2.9883e+06,1473.47 99550,2.98828e+06,1473.57 99550,2.98825e+06,1473.63 - To me it seems the easiest way to import these is the following: 1. v.in.ascii 2. Set the region to the extents of new vector file and the resolution to 25m. 3. v.to.rast. This 'works' but I get horizontal strips of no data in my raster DEM at 25m spacings. Also, when I look at the raster and the vector layer together, the vector points are not always on the edges of the raster cells. Can someone perhaps help me to fix this? g.region: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: 2988300 south: 2959850 west: 74300 east: 99550 nsres: 25 ewres: 25 rows: 1138 cols: 1010 cells: 1149380 v.info: ++ | Layer: dem_2628cc_...@c83 | Mapset: C83 | Location:sa_tm_19deg_E | Database:C:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: Administrator | Organization: | Source date: Thu May 13 16:06:07 2010 || | Type of Map: vector (level: 2) | | Number of points: 1151529 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 0 | | Projection: Transverse Mercator | N: 2988300S: 2959850 | E: 99550W: 74300 | B: 1429.79T:1740.2 | | Digitization threshold: 0 | Comments: | ++ r.info: ++ | Layer:dem_2628cc_...@c83 Date: Thu May 13 16:15:50 2010 | Mapset: C83Login of Creator: Administrator | Location: sa_tm_19deg_E | DataBase: C:\Hanlie\grassdata | Title:Categories ( dem_2628cc_25m ) | Timestamp: none || | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1138 | Columns: 1010 | Total Cells: 1149380 |Projection: Transverse Mercator |N:2988300S:2959850 Res:25 |E: 99550W: 74300 Res:25 | Range of data:min = 1429.85 max = 1739.41 | | Data Source: |Vector Map: dem_2628cc_...@c83 in mapset C83 |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=dem_2628cc_...@c83 output=dem_2628cc_25m use=z\ | type=point,line,area layer=1 value=1 rows=4096 | ++ Thanks Hanlie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: dems from coordinate lists
Did you adjust the region also? To do this properly, you need to first scan the file using r.in.xyz to get the extents. I'd also use the shell style output. r.in.xyz -sg input=2628cc.ORT.xyz output=dem_2628cc_25m_xyz Should return something like: n=1000 s=500 w=500 e=1000 t=### b=### (you can ignore the t b output) This is the extent of the points, which are going to be the center of the output pixels. Since the region is defined as pixel edge, you need to expand the region 1/2 a pixel in each direction to ensure the points are in the center of the pixel. At 25m spacing, your region should be set as: g.region n=1012.5 s=487.5 w=487.5 e=1012.5 res=25 - Jamie On Thu, May 13, 2010 at 9:02 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi Jamie, I tried r.in.xyz: (r.in.xyz input=2628cc.ORT.xyz output=dem_2628cc_25m_xyz method=mean type=FCELL x=1 y=2 z=3 zscale=1.0 percent=100) but I still get the rows of no data in the result. Hanlie You want the points to represent cell centers so you need to expand your region 1/2 pixel in all four directions. Also look at using r.in.xyz, it works directly on this type of data. Jamie On Thu, May 13, 2010 at 7:26 AM, Hanlie Pretorius hanlie.pretor...@gmail.com wrote: Hi, I've obtained DEMs in text files with columns X, Y and Z at 25m spacing. The first three entries in the text file are: - X,Y,Z 99550,2.9883e+06,1473.47 99550,2.98828e+06,1473.57 99550,2.98825e+06,1473.63 - To me it seems the easiest way to import these is the following: 1. v.in.ascii 2. Set the region to the extents of new vector file and the resolution to 25m. 3. v.to.rast. This 'works' but I get horizontal strips of no data in my raster DEM at 25m spacings. Also, when I look at the raster and the vector layer together, the vector points are not always on the edges of the raster cells. Can someone perhaps help me to fix this? g.region: projection: 99 (Transverse Mercator) zone: 0 datum: ** unknown (default: WGS84) ** ellipsoid: wgs84 north: 2988300 south: 2959850 west: 74300 east: 99550 nsres: 25 ewres: 25 rows: 1138 cols: 1010 cells: 1149380 v.info: ++ | Layer: dem_2628cc_...@c83 | Mapset: C83 | Location:sa_tm_19deg_E | Database:C:\Hanlie\grassdata | Title: | Map scale: 1:1 | Map format: native | Name of creator: Administrator | Organization: | Source date: Thu May 13 16:06:07 2010 || | Type of Map: vector (level: 2) | | Number of points: 1151529 Number of areas: 0 | Number of lines:0 Number of islands:0 | Number of boundaries: 0 Number of faces: 0 | Number of centroids:0 Number of kernels:0 | | Map is 3D: Yes | Number of dblinks: 0 | | Projection: Transverse Mercator | N: 2988300S: 2959850 | E: 99550W: 74300 | B: 1429.79T:1740.2 | | Digitization threshold: 0 | Comments: | ++ r.info: ++ | Layer:dem_2628cc_...@c83 Date: Thu May 13 16:15:50 2010 | Mapset: C83Login of Creator: Administrator | Location: sa_tm_19deg_E | DataBase: C:\Hanlie\grassdata | Title:Categories ( dem_2628cc_25m ) | Timestamp: none || | | | Type of Map: raster Number of Categories: 0 | Data Type:DCELL | Rows: 1138 | Columns: 1010 | Total Cells: 1149380 |Projection: Transverse Mercator |N:2988300S:2959850 Res:25 |E: 99550W: 74300 Res:25 | Range of data:min = 1429.85 max = 1739.41 | | Data Source: |Vector Map: dem_2628cc_...@c83 in mapset C83 |Original scale from vector map: 1:1 | | Data Description: |generated by v.to.rast | | Comments: |v.to.rast input=dem_2628cc_...@c83 output=dem_2628cc_25m use=z\ | type=point,line,area layer=1 value=1 rows=4096 | ++ Thanks Hanlie ___ grass-user mailing list
Re: [GRASS-user] can't ingest an E00 file using r.in.gdal, or anything else?
Hey Jeff, This is a valid file, though not well supported by software other than ESRI. E00 files most commonly hold vector formats, but can be used to contain raster data also. AFAIK, there is no way to import these directly into GRASS, it only supports the vector version with v.in.e00 For this reason, I keep a little ESRI utility around just in case: http://www.esri.com/apps/products/download/?downloadid=175 http://www.esri.com/apps/products/download/?downloadid=175It's a Windows binary, and requires registration, but it works. I've successfully run it using Wine on Linux, probably will work on a Mac also. I reads the e00 file and dumps an arcinfo binary grid where you specify. Hope this helps. Jamie On Sat, Oct 17, 2009 at 8:07 AM, Jeff Hamann jeff.ham...@forestinformatics.com wrote: I'm not sure why this is not working, and I'm having a little trouble getting a satisfactory answer on the web. I'm trying to import an e00 file (arc/info grid) using r.in.gdal. For example, GRASS 6.4.0svn (himom):~/Documents/rufus/incoming/wadnr r.in.gdal input=lewis_dem30.e00 output=lewisdem location=himom2 --verbose ERROR 4: `lewis_dem30.e00' not recognised as a supported file format. from outside GRASS, Jeff-Hamanns-MacBook-Pro:lewis_dir hamannj$ gdalinfo lewis_dem30.e00 ERROR 4: `lewis_dem30.e00' not recognised as a supported file format. gdalinfo failed - unable to open 'lewis_dem30.e00'. Jeff-Hamanns-MacBook-Pro:lewis_dir hamannj$ and the file looks like this: Jeff-Hamanns-MacBook-Pro:lewis_dir hamannj$ head lewis_dem30.e00 EXP 0 /AM/PROCESS/DATASERVICES/LEWIS_DIR/LEWIS_DEM30.E00 GRD 2 4860 1539 1-0.2147483647E+10 0.104355E+03 0.104355E+03 0.91924128086325E+06 0.38480838058312E+06 0.14264065808632E+07 0.54541072558312E+06 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 -2147483647 Jeff-Hamanns-MacBook-Pro:lewis_dir hamannj$ Any hints? What kind of file is this, really? Thanks, Jeff. Jeff Hamann, PhD PO Box 1421 Corvallis, Oregon 97339-1421 541-754-2457 jeff.hamann[at]forestinformatics[dot]com http://www.forestinformatics.com ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] nnbathy
Reviving this old thread - I just noticed the other day that the nn code has been checked into Google Code using an MIT license. Would it be possible to integrate directly now? I find r.surf.nnbathy to be one of the more useful interpolation routines. - Jamie On Tue, Sep 23, 2008 at 8:50 AM, Dylan Beaudette dylan.beaude...@gmail.comwrote: On Tue, Sep 23, 2008 at 6:12 AM, Moritz Lennert mlenn...@club.worldonline.be wrote: On 23/09/08 14:06, Moritz Lennert wrote: On 23/09/08 13:49, Moritz Lennert wrote: I have also put it on my machine, in case the above link doesn't work anymore: http://geog-pc40.ulb.ac.be/grass/nnbathy/ As a follow-up: I have sent a mail to Pavel Sakov to find out where the code is available now. Pavel sent me the code. It is available at above address. We will see whether he will set up a new web site, or whether we will integrate the code into the Add-Ons SVN. Moritz I have not seen the original code / license, but would it be possible to directly integrate the code into GRASS? This would remove the dependency of an external nnbathy library. Cheers, Dylan ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] mosaic raster files
If you have GDAL 1.6 installed, you can use gdalbuildvrt to create a VRT of the files, and then link to this using r.external. This way you have a mosaic without duplicating the data. On Thu, Jul 2, 2009 at 7:53 AM, Milton Cezar Ribeiro miltinho.astrona...@gmail.com wrote: Hi Guys, I will try first download the latest SVN. In fact I am running a compiled version of 6.4R4, and I remember, on that time, few programs (like r.patch if I am not wrong) not compiled properly. Case it not work, I try r.series. Thanks a lot, milton 2009/7/2 Hamish hamis...@yahoo.com Milton Cezar Ribeiro wrote: I just downloaded some tiles from GDEM/NASA and imported on GRASS. Now I would like to generate a mosaic with the tiles, but unfortunatelly r.patch ( http://grass.itc.it/gdp/html_grass63/r.patch.html ) is not compiled for GRASS 6.4R4 (windows / vista). Which installer? really there is nothing special about r.patch so it is a bad sign that it would be missing. what if you install RC5 or the SVN snapshot binary from June 2? So is there a way of I join a list of rasters on a unique mosaic file? r.series and r.mapcalc should both be able to do the job, but not as simply. use 'g.region rast=tile1,tile2,tile3,tile5' first to set the computational region to the outside bounds of all the tiles listed. Hamish ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] python commands on grass.
Try using help() from within the python shell.
python
import grass
help(grass)
-Jamie
On Tue, May 5, 2009 at 11:52 AM, Milton Cezar Ribeiro
miltinho.astrona...@gmail.com wrote:
Dear all,
is there a way of I know the commands available on python grass library?
Checking the codes available at
https://svn.osgeo.org/grass/grass/trunk/script
(thanks Martin) I saw the following commands:
grass.read_command
grass.parse_key_val
grass.run_command
grass.message
grass.fatal
By the way, how can I know the methods / parameters for each command?
Meanwhile, I can run g.list with the command bellow, but I can't
deal with the name of each of my raster map. In fact I receive the
same output as in gis.m console, but I would like only receive a list of
maps as output.
rasterfiles = grass.read_command('g.list',type='rast')
Cheers
milton
___
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
___
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] problem in projection of shape files
I would use ogr2ogr http://www.gdal.org/ogr2ogr.html from GDAL to reproject using a proj4 string. For example to wgs84: ogr2ogr -s_srs +proj=tmerc +lat_0=0 +lon_0=90 +k=0.9996 +x_0=50 +y_0=-200 +a=6377276.345 +b=6356075.41314024 +units=m +no_defs -t_srs wgs84 out.shp in.shp -Jamie On Mon, Feb 23, 2009 at 1:19 PM, Zahid Parvez grassgisbanglad...@gmail.comwrote: hi I have a shape file which I have to import in GRASS. Shape file is consists of a shape file, a dbf file and a prj file. The prj file has following information. PROJCS[BTM, GEOGCS[GCS_Everest_1830, DATUM[D_Everest_1830, SPHEROID[Everest_1830,6377299.36,300.8017]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.0174532925199433]], PROJECTION[Transverse_Mercator], PARAMETER[False_Easting,50.0], PARAMETER[False_Northing,-200.0], PARAMETER[Central_Meridian,90.0], PARAMETER[Scale_Factor,0.9996], PARAMETER[Latitude_Of_Origin,0.0], UNIT[Meter,1.0]] I will be very grateful if anyone could help me to find the ESPG code and any other way to re-project the shape file and import into GRASS. Best Md. Zahid Parvez ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] saving patch areas on raster maps
You can use r.recode to do this, but you'll need to reformat the output from r.stats. Example: *r.stats -a -n raster* 0 1234 1 2345 2 3456 .. This needs to be: 0:0:1234:1234 1:1:2345:2345 2:2:3456:3456 .. Save this to a txt file and use it with the rules flag of r.recode. If you have a large amount of categories, I would suggest using awk or perl to reformat the r.stats output. -Jamie On Fri, Nov 14, 2008 at 8:22 AM, Milton Cezar Ribeiro [EMAIL PROTECTED] wrote: Dear all, I generated a clump map from my cover map, and computed the area of each patch using r.stats. Now I need generate a map which values is the area for each clump. Suggestions are welcome. Best regards, miltinho astronatua brazil ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Re: r.resamp.stats max method not working as expected
I think 7 is fine. My current approach is working well enough. Thanks! On Thu, Oct 16, 2008 at 2:34 AM, Markus Neteler [EMAIL PROTECTED] wrote: On Wed, Oct 15, 2008 at 5:39 PM, Glynn Clements [EMAIL PROTECTED] wrote: Glynn Clements wrote: ... Implemented in r33888 (in 7.0). The aforementioned comment in the Wiki still stands. min/max aren't actually weighted, so the behaviour of those aggregates with a perfectly aligned grid is unstable (i.e. the exact result may depend upon miniscule rounding errors) even when -w is used. Jamie, do you want need this in 6.4 (hence, backport)? Or are you fine with GRASS 7? According to Glynn it's a slightly incompatible change which isn't 100% in line with our backport policy. Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] Re: r.resamp.stats max method not working as expected
Okay, I did find verification that this doesn't work in r.resamp.stats: *[1] The min and max aggregates can't use weights, so -w has no effect for those.* From http://grass.osgeo.org/wiki/GRASS_raster_semantics - a good read for anyone working with raster data. -Jamie On Mon, Oct 13, 2008 at 12:54 PM, Jamie Adams [EMAIL PROTECTED] wrote: Hello all, When I use r.resamp.stats to resample a grid of 0s 1s using the maximum option and the -w flag, it doesn't seem to consider the edge pixels that overlap. It returns the same result as without the -w flag. I've gotten around this by using the default method and setting anything 0 to 1. Any reason for this behavior, other than being a bug? -Jamie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Re: r.resamp.stats max method not working as expected
Very cool. I went ahead and filed a enhancement ticket for this in trac: http://trac.osgeo.org/grass/ticket/337 On Tue, Oct 14, 2008 at 2:34 AM, Glynn Clements [EMAIL PROTECTED]wrote: Jamie Adams wrote: Okay, I did find verification that this doesn't work in r.resamp.stats: *[1] The min and max aggregates can't use weights, so -w has no effect for those.* Hmm. There's no fundamental reason why I can't add weighted versions of the min/max aggregates. Obviously, these would ignore the weights, but it would allow r.resamp.stats to include the source cell in the calculation for each destination cell which overlaps it. -- Glynn Clements [EMAIL PROTECTED] ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.in.gdal global raster incorrect extents
but why then does it work fine for me? I downloaded a 180W tile, that looks fine too: ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM3/Australia/S20W180.hgt.zip Origin = (-180.0004166,-18.9995834) Pixel Size = (0.0008333,-0.0008333) I grabbed the same tile and imported it into GRASS using r.in.srtm and am having the same issue: gdalinfo S20W180.hgt Origin = (-180.0004166,-18.9995834) Pixel Size = (0.0008333,-0.0008333) r.in.srtm input=S20W180 output=S20W180 cat /location/mapset/cellhd/S20W180 proj: 3 zone: 0 north: 18:59:58.5S south: 20:00:01.5S east: 178:59:58.5W west: 179:59:58.5E cols: 1201 rows: 1201 e-w resol: 0:00:03 n-s resol: 0:00:03 format: 3 compressed: 1 gdalinfo /location/mapset/cellhd/S20W180 Origin = (179.9995834,-18.9995834) Pixel Size = (0.0008333,-0.0008333) r.out.gdal outputs a GeoTiff consistent with this gdalinfo output - so shifted 180E. where do your .tif files come from? I wrote a python script that creates these as equal divisions of an input area (in this case the world). I wanted to use them along with srtm, so I used the same res and pixel registration scheme. I looked into it a bit further and noticed that my script is using a truncated version of the res variable (initialized to 1/1200): geotransform: (-180.000416667, 0.00083329, 0.0, 10.000417, 0.0, -0.00083329) cat cellhd/W180N10 proj: 3 zone: 0 north: 10:00:01.5N south: 10:00:01.5S east: 159:59:58.51W west: 179:59:58.49E cols: 24001 rows: 24001 e-w resol: 0:00:03 n-s resol: 0:00:03 format: 0 compressed: 1 I fixed this and the precision issue went away. But as with the srtm tile, it doesn't fix the e-w region issue. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] r.in.gdal global raster incorrect extents
Hello all, When I import a global elevation dataset with extents 180w 90n 180e 90s using r.in.gdal, the GRASS raster shows up with w-e extents both showing 180E. While this doesn't seem to be affecting raster calculations, afaik, it's causing issues with QGIS. It displays it with it's origin at 180E and extends to 540E. Has anyone run into this problem? *gdalinfo:* Driver: GTiff/GeoTIFF Files: 2_00-res/out-2_00-nozero-res.tif Size is 43200, 21600 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.2572235630016, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]] Origin = (-180.0124970,90.000) Pixel Size = (0.008,-0.008) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left (-180.000, 90.000) (180d 0'0.00W, 90d 0'0.00N) Lower Left (-180.000, -90.000) (180d 0'0.00W, 90d 0'0.00S) Upper Right ( 180.000, 90.000) (180d 0'0.00E, 90d 0'0.00N) Lower Right ( 180.000, -90.000) (180d 0'0.00E, 90d 0'0.00S) Center ( -0.000, 0.000) ( 0d 0'0.00W, 0d 0'0.00N) Band 1 Block=43200x1 Type=Float32, ColorInterp=Gray NoData Value=10 * r.info:* ++ | Layer:2_00_res Date: Fri Oct 10 17:45:37 2008 | Mapset: PERMANENT Login of Creator: jamie | Location: Bathy | DataBase: /work/grass | Title: ( 2_00_res ) | Timestamp: none | | | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 21600 | Columns: 43200 | Total Cells: 93312 |Projection: Latitude-Longitude |N:90NS:90S Res: 0:00:30 |E: 180EW: 180E Res: 0:00:30 | Range of data:min = -10814.00 max = -0.01 | | Data Description: |generated by r.in.gdal | | Comments: |r.in.gdal input=out-2_00-nozero-res.tif output=2_00_res | ++ ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.in.gdal global raster incorrect extents
Thanks for the suggestions. I opened the image up in python using gdal, looked at the geotransform and found that the image is ever-so-slightly shifted west (-180.012497). Hence the negative W-E center. Maybe this is why GRASS thinks it has a 180E origin also, though I don't know. I checked the region - same problem - then set it manually and re-imported. Same results. In the end I was able to assign correct extents using r.region. Though the boundary isn't exactly at 180W, I think nine 0s after the decimal point is a bit much. :-) Thanks again. -Jamie On Mon, Oct 13, 2008 at 1:35 PM, Nikos Alexandris [EMAIL PROTECTED] wrote: On Mon, 2008-10-13 at 12:41 -0700, Jamie Adams wrote: Hello all, When I import a global elevation dataset with extents 180w 90n 180e 90s using r.in.gdal, the GRASS raster shows up with w-e extents both showing 180E. While this doesn't seem to be affecting raster calculations, afaik, it's causing issues with QGIS. It displays it with it's origin at 180E and extends to 540E. Has anyone run into this problem? gdalinfo: Driver: GTiff/GeoTIFF Files: 2_00-res/out-2_00-nozero-res.tif Size is 43200, 21600 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.2572235630016, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]] Origin = (-180.0124970,90.000) Pixel Size = (0.008,-0.008) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left (-180.000, 90.000) (180d 0'0.00W, 90d 0'0.00N) Lower Left (-180.000, -90.000) (180d 0'0.00W, 90d 0'0.00S) Upper Right ( 180.000, 90.000) (180d 0'0.00E, 90d 0'0.00N) Lower Right ( 180.000, -90.000) (180d 0'0.00E, 90d 0'0.00S) Center ( -0.000, 0.000) ( 0d 0'0.00W, 0d 0'0.00N) Band 1 Block=43200x1 Type=Float32, ColorInterp=Gray NoData Value=10 r.info: ++ | Layer:2_00_res Date: Fri Oct 10 17:45:37 2008 | Mapset: PERMANENT Login of Creator: jamie | Location: Bathy | DataBase: /work/grass | Title: ( 2_00_res ) | Timestamp: none | | | Type of Map: raster Number of Categories: 255 | Data Type:FCELL | Rows: 21600 | Columns: 43200 | Total Cells: 93312 |Projection: Latitude-Longitude |N:90NS:90S Res: 0:00:30 |E: 180EW: 180E Res: 0:00:30 | Range of data:min = -10814.00 max = -0.01 | | Data Description: |generated by r.in.gdal | | Comments: |r.in.gdal input=out-2_00-nozero-res.tif output=2_00_res | ++ 1. I am not 100% sure but I think the line Center ( -0.000, 0.000) ( 0d 0'0.00W, 0d 0'0.00N) in the geo-metadata should read enter ( 0.000, -0.000) ( 0d 0'0.00W, 0d 0'0.00N) (Could it be that the file is not ok?) 2. What doeas g.region -p say? If you get the same value for east and west it's wrong I guess. west should be, well... west: 180W and east: 180E. 3. Try to set g.region e=180 s=-90 w=-180 n=90 -p. Re-import again and check for errors. Regards, Nikos ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.in.gdal global raster incorrect extents
Sure, python w/ gdal is great for doing these sanity checks. Assuming you
are using a FWTools install or have the python-gdal plugin:
python
import gdal
or newer version:
from osgeo import gdal
image = gdal.Open('/the/path/to/my/raster.tif',gdal.GA_ReadOnly)
or if your version of gdal has Grass support
image =
gdal.Open('/grassdir/location/mapset/cellhd/rastername',gdal.GA_ReadOnly)
image.GetGeoTransform()
(-180.0, 0.0083332, 0.0, 90.0, 0.0, -0.0083332)
In my case, the source GeoTiff had a geotransform of:
(-180.012497, 0.0083297, 0.0, 90.0, 0.0,
-0.0083297)
One strange thing was, the original grass import (west coordinate of 180E)
had the same geotransform as the corrected raster I ran r.region on.
According to gdal, they were the same.
Looks like I need to file a bug report.
Enjoy.
-Jamie
On Mon, Oct 13, 2008 at 4:18 PM, Nikos Alexandris
[EMAIL PROTECTED] wrote:
On Mon, 2008-10-13 at 16:13 -0700, Jamie Adams wrote:
Thanks for the suggestions. I opened the image up in python using
gdal, looked at the geotransform and found that the image is
ever-so-slightly shifted west (-180.012497). Hence the
negative W-E center. Maybe this is why GRASS thinks it has a 180E
origin also, though I don't know.
I checked the region - same problem - then set it manually and
re-imported. Same results. In the end I was able to assign correct
extents using r.region. Though the boundary isn't exactly at 180W, I
think nine 0s after the decimal point is a bit much. :-)
Thanks again.
-Jamie
Jamie,
would you share some python command examples (using gdal)? I am learning
python stuff (although very slowly) and I would like to have some
examples as the one you did.
Thank you, Nikos
___
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] r.in.gdal global raster incorrect extents
On Mon, Oct 13, 2008 at 6:08 PM, Hamish [EMAIL PROTECTED] wrote: Jamie Adams wrote: In my case, the source GeoTiff had a geotransform of: (-180.012497, 0.0083297, 0.0, 90.0, 0.0, -0.0083297) One strange thing was, the original grass import (west coordinate of 180E) had the same geotransform as the corrected raster I ran r.region on. It may have on the surface (ie what the program output reports: up to the programmer defined number of decimal places), but internally the stored value may be ever so slightly different. the actual value used for grass rasters is in $MAPSET/cellhd/$MAPNAME it is still not ok after r.region? My global 0.5 min grid is working, but I've discovered problems with other data sets fyi, you can fix/reset the broken GeoTiff with gdal_translate -a_ullr Yes, I'm aware of that - just questioning what a broken file is (see below) According to gdal, they were the same. Using python-gdal ie gdalinfo r.info? or python-gdal? Looks like I need to file a bug report. with QGIS? what is the problem besides a slightly-broken input file? I'm not sure why this is considered a broken input file. Yes it has a slight shift to the west, but all of the cell centers lie between 180w 90n 180e 90s. Yes, this amount of shift is trivial and can be easily fixed, but it could be considerably worse. Take for example some other raster data I have that is registered up to 180w, but has the pixel centers aligned to the grid (i think srtm is this way). I do a r.in.gdal, g.region rast=in.tif, r.out.gdal and end up with a completely different registration. gdalinfo W180N10.tif Driver: GTiff/GeoTIFF Files: W180N10.tif Size is 24001, 24001 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.2572235630016, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]] Origin = (-180.0004196,10.0004170) Pixel Size = (0.0008333,-0.0008333) Metadata: AREA_OR_POINT=Area Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND Corner Coordinates: Upper Left (-180.0004167, 10.0004167) (180d 0'1.50W, 10d 0'1.50N) Lower Left (-180.0004167, -10.0004167) (180d 0'1.50W, 10d 0'1.50S) Upper Right (-159.9995833, 10.0004167) (159d59'58.50W, 10d 0'1.50N) Lower Right (-159.9995833, -10.0004167) (159d59'58.50W, 10d 0'1.50S) Center (-170.000, 0.000) (170d 0'0.00W, 0d 0'0.00N) Band 1 Block=24001x1 Type=Byte, ColorInterp=Gray r.in.gdal input=W180N10.tif output=W180N10 r.info -g W180N10 north=10:00:01.5N south=10:00:01.5S east=159:59:58.51W west=179:59:58.49E g.region rast=W180N10 r.out.gdal input=W180N10 output=W180N10_v2.tif gdalinfo W180N10_v2.tif Driver: GTiff/GeoTIFF Files: gdalout.tif Size is 24001, 24001 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.2572235630016, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]] Origin = (179.99958055540,10.0004166) Pixel Size = (0.0008333,-0.0008333) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 179.9995833, 10.0004167) (179d59'58.50E, 10d 0'1.50N) Lower Left ( 179.9995833, -10.0004167) (179d59'58.50E, 10d 0'1.50S) Upper Right ( 200.000, 10.000) (200d 0'1.50E, 10d 0'1.50N) Lower Right ( 200.000, -10.000) (200d 0'1.50E, 10d 0'1.50S) Center ( 190.000, 0.000) (190d 0'0.00E, 0d 0'0.01N) That seems like a GRASS, maybe GDAL issue. -Jamie Hamish ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] get modified date for raster layer
Hello all, I've got over 14k rasters in several mapsets that I need to store modified date in a db for tracking purposes. I've got the storage and retrieval worked out, but I don't see an clean (easy) way to get the modified date via a python script. So far I've come up with: - parse it out of r.info - get the modified date in the cellhd folder using os.stat The problems with these approaches: r.info - requires parsing, and I need to reformat the date in a way that's easily comparable. I'd prefer not to rely on parsing anyway. Would be nice to have shell script output option. os.stat (python) - has a nice date integer format, but it requires that I have the full path to the cellhd folder. Database and location are easy via g.gisenv, but getting the mapset on an individual layer basis isn't completely straightforward. I can use use g.mlist -m rast pattern=filename, but that essentially greps it out of a full list. Is there a more efficient way of getting the modified date of a raster layer? If not, getting the full path instead? Thanks, Jamie Adams ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] get modified date for raster layer
On Tue, Sep 2, 2008 at 6:42 PM, Hamish [EMAIL PROTECTED] wrote: Jamie Adams [EMAIL PROTECTED] I've got over 14k rasters in several mapsets that I need to store modified date in a db for tracking purposes. I've got the storage and retrieval worked out, but I don't see an clean (easy) way to get the modified date via a python script. Is there a more efficient way of getting the modified date of a raster layer? If not, getting the full path instead? MAP=roads eval `g.findfile element=hist file=$MAP` head -n 1 $file It's stored as a text string; what you see is what you get. see also r.timestamp, but that is for map metadata not file creation timestamp. os.stat (python) - has a nice date integer format, but it requires that I have the full path to the cellhd folder. use g.findfile for that. Ah okay, I never really understood that command. g.findfile element=cell file=filename | grep ^file returns what I need. Database and location are easy via g.gisenv, but getting the mapset on an individual layer basis isn't completely straightforward. I can use use g.mlist -m rast pattern=filename, but that essentially greps it out of a full list. I don't fully understand what you mean by those two. Yeah, it's a bit tough to explain with my processing setup. It's a non-issue now, g.findfile solved it. Thanks! Hamish ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Re: make a raster index in GRASS
Thanks for the tip. I ended up using FWTools on linux w/ Grass support and gdaltindex worked fine. -Jamie On Wed, Jul 16, 2008 at 6:07 AM, G. Allegri [EMAIL PROTECTED] wrote: In case you have compiled gdal with grass support you could access grass rasters directly from gdaltindex, pointing to the cellhd files. 2008/7/16 Ivan Shmakov [EMAIL PROTECTED]: Jamie Adams [EMAIL PROTECTED] writes: Hello all, I'd like to generate a raster index polygon file, just like gdaltindex does, but using GRASS rasters as input. Is there a command I'm overlooking? It may depend on the goal, but v.in.region(1) may be of some use. If not, any ideas on how to do this? Iterating over the set of rasters can be implemented using g.mlist(1) and the standard Shell `while' and `read' commands, like: $ g.mlist type=rast pattern=2008-\*-foo \ | while read r ; do \ : ... do something... ; \ done And the vectors can be concatenated using v.patch(1), like: v.in.region output=tmp_vector v.patch -a output=resulting_vector input=tmp_vector ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] make a raster index in GRASS
Hello all, I'd like to generate a raster index polygon file, just like gdaltindex does, but using GRASS rasters as input. Is there a command I'm overlooking? If not, any ideas on how to do this? -Jamie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] SRTM and r.fillnulls - using neighboring tiles
Hello all, I'm trying to devise an efficient way to fill voids in SRTM, using neighboring tiles when there is edge nodata. In the end, I want to have tiles matching the input, just lacking any voids. My initial thought is to write a simple gdal python script to determine which tiles have edge nodata and only process those. I'll then write a script to: - determine which tiles are neighbors - merge those into a single temp dem - set the region to the desired tile - r.fillnull the temp file, creating an output with the correct extents This should work, though I'm thinking there must be a more efficient way. Anyone have suggestions? Thanks! Jamie ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Orthorectification of Historic Aerial Photographs
Would it be possible to calculate these parameters with a panoramic creation package like panotools? http://wiki.panotools.org/Lens_correction_model It calculates a value for field of view, but I didn't notice an actual focal length parameter. Maybe this is easy to associate? One benefit is that you are calculating using image overlap, which may be easier than trying to reference to a recent photo. -Jamie On Thu, Apr 10, 2008 at 9:53 AM, Vincent Bain [EMAIL PROTECTED] wrote: May it be a stupid suggestion, but did you try to contact anyone at the USGS ? perhaps they maintain data from these old times... (here in France IGN is always helpful, and can provide answers to such historical requests) Vincent. Le jeudi 10 avril 2008 à 18:33 +0200, Markus Neteler a écrit : On Thu, Apr 10, 2008 at 5:29 PM, Shaun Melissa Busler [EMAIL PROTECTED] wrote: Unfortunately, the focal length of the camera is not on the photo border. The photo must be too old. Do you have any recommendations on what focal length I should try? It was taken in 1939 for the USDA Agricultural Adjustment Administration. You could try 150mm or 305mm. In fact, i.ortho.photo is quite robust, just try and hope that it converges. (and/or: Maybe someone on the PROJ4 list knows.) Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
