[GRASS-dev] how to sample a series at one location?

[Newcomb, Doug]

Michael,
You could use v.what.rast  in a python script , iterating the raster layers
with an sqlite database back end.  I think that would go up to 2000 columns
for each point you have.

Alternatively, you could use a bit of python with gdal.  I was trying to do
something similar in GRASS, to change the Z value of each point in a text
LiDAR file, from absolute above sea level to relative the elevation of the
ground.  For 25.5 billion text points and a Statewide 20 ft elevation grid
(6.5 ? billion cells) , it was a bit slow using r.what.  So I converted the
LiDAR data to (7)  3.3 billion point  LAS format files and exported the
GRASS layer to an Erdas imagine format file and wrote the following ugly
python script:

#!/usr/bin/python
import os,string,glob,re,gdal
from liblas import file
from liblas import header
from liblas import point
from gdalconst import *
h=header.Header()
#enter the LAS point file name
infile=raw_input("Enter the input lidar data points file: ")
#Hardcoded edras imagine file, you will have to use an array for the
different data layer names
imgfile="/gisdata2/raster/allnc_20ft_el.img"
#print "suggest /gisdata2/raster for output dir\n"
inarr=infile.split('.')
#This sets the LAS output file, substitute your output text file instead
outfil=inarr[0]+"_norm.las"
#outfil=raw_input("Enter output text file name: ")
#This part reads the LAS file, if you
l=file.File(infile,mode='r')
#Outputs the LAS file
lout=file.File(outfil,mode='w',header=h)
# register all of the drivers, hopefully your gdal speaks GRASS
gdal.AllRegister()
#opening and closing the image layers might take some time if you are
reading thousands of images
ds=gdal.Open(imgfile,GA_ReadOnly)
if ds is None:
    print 'Could not open image'
    sys.exit(1)
# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
# get georeference info, not sure how this would work for GRASS data layers
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transformAsArray(xOffset, yOffset, 1, 1)
pixelWidth = transform[1]
pixelHeight = transform[5]
for p in l:
    x=float(p.x)
    y=float(p.y)
    z=float(p.z)
    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    band = ds.GetRasterBand(1) # 1-based index 0? 1?
    data = band.Readr(value) :continue
    value = data[0,0]
    #print value,"11","\n"
    if "nan" in st[3]
    znorm = z-value
    #print znorm,"\n"
    pt=point.Point()
    pt.x=p.x
    pt.y=p.y
    pt.z=znorm
    lout.write(pt)

l.close()
lout.close()
#25561312019 points in allreturns

This managed to process everything over a weekend( in 7 parallel threads),
which was fast enough for me at the time.  Approaching your problem , if
your data layers are all n GRASS and your version of gdal can read GRASS
data layers,  I would grab the list of GRASS layers via glob and iterate
the layers , writing the name of the raster layer and the value of the
raster layer at the point coordinates out to a text file as you state
below.

Hope this helps,
Doug


On Fri, Feb 22, 2013 at 1:23 PM, Michael Barton <Michael.Barton at asu.edu>wrote:

> But I want to do it with a time series of hundreds or thousands of maps.
>
> Michael
> ______________________________
> C. Michael Barton
> Director, Center for Social Dynamics & Complexity
> Professor of Anthropology, School of Human Evolution & Social Change
> Arizona State University
> Tempe, AZ  85287-2402
> USA
>
> voice:  480-965-6262 (SHESC), 480-727-9746 (CSDC)
> fax:          480-965-7671(SHESC), 480-727-0709 (CSDC)
> www:    http://csdc.asu.edu, http://shesc.asu.edu
>                 http://www.public.asu.edu/~cmbarton
>
> On Feb 22, 2013, at 10:53 AM, Markus Neteler <neteler at osgeo.org>
>  wrote:
>
> > On Fri, Feb 22, 2013 at 6:41 PM, Michael Barton <Michael.Barton at asu.edu>
> wrote:
> >> Is there tool somewhere, including in the new temporal GIS modules, to
> >> sample the value of a raster series at one cell location? I'd like to
> get
> >> text something like this for a cell that I specify with xy coordinates
> or a
> >> cat:
> >>
> >> map1, value 1
> >> map2, value 2
> >> map3, value 3
> >> map4, value 4
> >> map5, value 5
> >> …
> >>
> >
> > I do such queries with r.what which accepts multiple input:
> >
> > http://grass.osgeo.org/grass64/manuals/r.what.html
> >
> > Markus
>
> _______________________________________________
> grass-dev mailing list
> grass-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-dev
>



-- 
Doug Newcomb
USFWS
Raleigh, NC
919-856-4520 ext. 14 doug_newcomb at fws.gov
---------------------------------------------------------------------------------------------------------
The opinions I express are my own and are not representative of the
official policy of the U.S.Fish and Wildlife Service or Dept. of the
Interior.   Life is too short for undocumented, proprietary data formats.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20130222/0e5ceea9/attachment-0001.html>