[GRASS-user] making fcell raster maps using climate data

GRASS GRASS user
1

Dear List,

I have an Excel file that contains precipitation data in 13 columns (mean values for each month and annual average) from a region between certain years (in rows): eventually showing temporal changes in precipitation values. The data are taken at a weather station whose coordinates are known. I want to use this file in order to make an fcell raster map in GRASS. After some research I found few leads and mostly on interpolation of data. However, I am confused about the actual process turning the data in Excel into a raster map. I will greatly appreciate any suggestions as to which module to use and how to arrange the data.

Thank you,

BÜLENT

2

Hi Bulent,

My first idea - but I may be biased as I come from the R world - would
be to make the imports steps of your data from R. R has a XLS
reader(but I usually find it safer to first export your spreadsheet to
a CSV file).

Hope this helps,

Pierre

2011/5/30 Bulent Arikan <bulent.arikan@gmail.com>:

Dear List,
I have an Excel file that contains precipitation data in 13 columns (mean
values for each month and annual average) from a region between certain
years (in rows): eventually showing temporal changes in precipitation
values. The data are taken at a weather station whose coordinates are known.
I want to use this file in order to make an fcell raster map in GRASS. After
some research I found few leads and mostly on interpolation of data.
However, I am confused about the actual process turning the data in Excel
into a raster map. I will greatly appreciate any suggestions as to which
module to use and how to arrange the data.

Thank you,
--
BÜLENT

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
Scientist
Landcare Research, New Zealand

3

Bulent,

Others will probably have different suggestions, but first you must get your data into GRASS. Both involve first using v.in.ascii to import (assuming you have latitude-longitude coordinates) the identifiers and locations (lat-long) into GRASS. Then you could, say, bring the climate data into a sqlite database and then connect your imported IDs and locations you previously imported into GRASS to the sqlite database.

Alternatinely, and probably easier is to use v.in.ascii with the station IDs, lat & long values and the 13 climate values for each station as attributes to the vector points. See the v.in.ascii documentation for examples.

The next step is to do the raster spatial interpolation from the vector point data — there are many ways to do this within GRASS, using something as simple as inverse distance weighting to splines with tension or to use R and gstat within GRASS and use kriging.

Regards,
Tom

----- Original Message -----
From: Bulent Arikan <bulent.arikan@gmail.com>
Date: Sunday, May 29, 2011 7:58 pm
Subject: [GRASS-user] making fcell raster maps using climate data
To: grass-user@lists.osgeo.org

Dear List,

I have an Excel file that contains precipitation data in 13 columns (mean
values for each month and annual average) from a region between certain
years (in rows): eventually showing temporal changes in precipitation
values. The data are taken at a weather station whose coordinates are
known.
I want to use this file in order to make an fcell raster map in GRASS.
After
some research I found few leads and mostly on interpolation of data.
However, I am confused about the actual process turning the data in Excel
into a raster map. I will greatly appreciate any suggestions as to which
module to use and how to arrange the data.

Thank you,
--
BÜLENT
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org

4

Bulent,

Thank you for taking time and helping out. I attached a short version of
data in (.xls) format. The file contains climate model results between 0 and
800 BP, for every century, at a single weather station (Van-Turkey).

Thanks. This is much clearer :slight_smile:

When you say "reformatting" do you actually mean organizing data for importing as
a vector point?

Yes. Consider the following R session. I just exported your XLS data
sheet into CSV in LibreOffice.

df <- read.csv('Van_Precipitation.csv')
df

   UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May
1 502500 4398240 1661 0 35.35559 34.10211 41.08517 54.96244 53.24204
2 502500 4398240 1661 -100 35.35434 34.10950 41.11511 54.98155 51.01671
3 502500 4398240 1661 -200 34.95282 34.74009 43.04097 56.12360 45.02433
4 502500 4398240 1661 -300 35.35213 34.12072 41.16004 55.01019 48.52734
5 502500 4398240 1661 -400 35.20788 34.33214 41.78447 55.40257 18.66510
6 502500 4398240 1661 -500 35.10954 34.49606 42.30229 55.70860 49.07919
7 502500 4398240 1661 -600 35.06167 34.58229 42.58510 55.86880 48.44339
8 502500 4398240 1661 -700 35.00383 34.69651 42.97523 56.08147 45.40806
9 502500 4398240 1661 -800 35.01142 34.68574 42.94451 56.06451 45.70229
        Jun Jul Aug Sep Oct Nov Dec Annual
1 17.729643 22.593137 0.7567056 11.717915 54.03994 41.26782 31.98477 398.8373
2 18.018884 22.797248 2.0428841 11.654079 53.91107 41.37436 31.99786 398.3736
3 10.662323 9.111732 4.5243246 4.254332 44.57887 48.04743 33.02368 368.0845
4 18.508462 23.089403 10.5560756 11.613746 53.75689 41.53623 32.01759 405.2488
5 21.111523 21.222403 19.7168880 10.175261 50.41907 43.68269 32.32056 384.0406
6 19.838364 3.158531 0.5505311 10.543923 47.83534 45.58733 32.58478 376.7945
7 15.610073 17.096418 8.2630408 12.148564 46.61419 46.60779 32.73860 395.6199
8 9.915775 8.735210 5.0867794 1.947931 45.02865 47.92435 32.96591 365.7697
9 10.045168 9.128862 5.3451146 2.768823 45.19434 47.83308 32.94591 367.6698

Then I'll use the reshape 2 to change the structure of df:

library(reshape2)
res <- melt(df, id.vars=c(1:4), variable.name='month', value.name='model_output')
head(res, 15)

    UTM_E UTM_N ELEV DATE month model_output
1 502500 4398240 1661 0 Jan 35.35559
2 502500 4398240 1661 -100 Jan 35.35434
3 502500 4398240 1661 -200 Jan 34.95282
4 502500 4398240 1661 -300 Jan 35.35213
5 502500 4398240 1661 -400 Jan 35.20788
6 502500 4398240 1661 -500 Jan 35.10954
7 502500 4398240 1661 -600 Jan 35.06167
8 502500 4398240 1661 -700 Jan 35.00383
9 502500 4398240 1661 -800 Jan 35.01142
10 502500 4398240 1661 0 Feb 34.10211
11 502500 4398240 1661 -100 Feb 34.10950
12 502500 4398240 1661 -200 Feb 34.74009
13 502500 4398240 1661 -300 Feb 34.12072
14 502500 4398240 1661 -400 Feb 34.33214
15 502500 4398240 1661 -500 Feb 34.49606

Here, I got the value of the result model for each month in each year,
for each station.

The following step is interpolation. This is not clear what you want
to do, but I guess what you want is to generate one raster layer per
month and per year.

I propose to use the plyr library, which allows you to split your data
in function of one or more variables (for us, month and year), and
apply a function to it.

Disclaimer: the following script is untested (that's just the general
idea) and should be launched from a R session WITHIN GRASS. Also, I
did not include the projection of the variables "grid" and "data".

# Generate an interpolation grid (the same for each raster to be generated)
resolution <- 5 # your resolution
grid <- expand.grid(x = seq(min(df$UTM_E), max(df$UTM_E), by =
resolution), y = seq(min(df$UTM_N), max(df$UTM_N), by = resolution))
grid <- points2grid(grid)

# This is a function that generate a map from each "slice" of data
make_map <- function(data, grid){
  # dependencies
  require(sp)
  require(rgdal)
  require(spgrass6)
  require(gstat) # for kriging
  require(automap) # to automate kriging

  # converting the data into a SpatialPointsDataFrame
  coordinates(data) <- ~UTM_E + UTM_N

  # INSERT YOUR FAVORITE INTERPOLATION METHOD HERE
  res <- autoKrige(model_output ~ 1, data, grid)

  # write it as a raster here
  layer_name <- paste("model_output_", unique(data$month), "-",
unique(data$DATE), sep ="") # a name for the current layer
  writeRAST6(res, layer_name, zcol="var1.pred", ignore.stderr=TRUE,
useGDAL=TRUE)

  # return the result
  return(res)
}

# Here is the plyr call - split the data, apply the function
"make_map", which interpolates the data, creates a layer name and
writes it in GRASS. Finally, it writes back everything in "rasters",
which is a list.
library(plyr)
rasters <- dlply(.data = res, .variables = .(month, DATE), .fun =
make_map, grid = grid)

Again, untested and probably buggy. Probably others can do better and
quicker with a bit of awk magic.

Hope this helps,

Pierre

Cheers,
Bulent

On Sun, May 29, 2011 at 8:19 PM, Pierre Roudier <pierre.roudier@gmail.com>
wrote:

No you probably don't need to do so. What R allows you to do, is to
reformat your data frame so you do not need 400 files!

Could you post a (short) overview of the structure of your data? From
what I understand, I'd use R and the plyr package to reformat your
data before attempting an import.

Cheers,

Pierre

2011/5/30 Bulent Arikan <bulent.arikan@gmail.com>:
> Hi Pierre,
> Thank you for answering. I use R too, so your reply was definitely
> helpful.
> I am a bit confused –and, I do not mean that I am asking for an answer
> from
> you when writing this– about the steps involving GRASS. So far I tried
> "r.in.ascii" and "r.in.bin" without success.Simply, the CSV file I have
> has
> 400 rows (each row represents one year) and I am not sure if I need to
> make
> individual CSV files for each year (means 400 files!!!) before importing
> into GRASS.
> Thank you for writing.
> Bulent
> On Sun, May 29, 2011 at 5:40 PM, Pierre Roudier
> <pierre.roudier@gmail.com>
> wrote:
>>
>> Hi Bulent,
>>
>> My first idea - but I may be biased as I come from the R world - would
>> be to make the imports steps of your data from R. R has a XLS
>> reader(but I usually find it safer to first export your spreadsheet to
>> a CSV file).
>>
>> Hope this helps,
>>
>> Pierre
>>
>> 2011/5/30 Bulent Arikan <bulent.arikan@gmail.com>:
>> > Dear List,
>> > I have an Excel file that contains precipitation data in 13 columns
>> > (mean
>> > values for each month and annual average) from a region between
>> > certain
>> > years (in rows): eventually showing temporal changes in precipitation
>> > values. The data are taken at a weather station whose coordinates are
>> > known.
>> > I want to use this file in order to make an fcell raster map in
>> > GRASS.
>> > After
>> > some research I found few leads and mostly on interpolation of data.
>> > However, I am confused about the actual process turning the data in
>> > Excel
>> > into a raster map. I will greatly appreciate any suggestions as to
>> > which
>> > module to use and how to arrange the data.
>> >
>> > Thank you,
>> > --
>> > BÜLENT
>> >
>> > _______________________________________________
>> > grass-user mailing list
>> > grass-user@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/grass-user
>> >
>> >
>>
>>
>>
>> --
>> Scientist
>> Landcare Research, New Zealand
>
>
>
> --
> BÜLENT ARIKAN, PhD
> Postdoctoral Scholar
> Center for Social Dynamics and Complexity &
> School of Human Evolution and Social Change
> Arizona State University
> Tempe - AZ
> 85287-2402
>

--
Scientist
Landcare Research, New Zealand

--
BÜLENT ARIKAN, PhD
Postdoctoral Scholar
Center for Social Dynamics and Complexity &
School of Human Evolution and Social Change
Arizona State University
Tempe - AZ
85287-2402

--
Scientist
Landcare Research, New Zealand

5

Pierre,
I cannot thank you enough for all these detailed, step-by-step explanations.

No worries :wink:

I am familiar with "v.surf.rst" as an
interpolation method, so I will use that.

Then you can probably insert a system call to GRASS instead of my
kriging routine in make_maps().

Pierre

Thank you again,
Bulent

On Sun, May 29, 2011 at 11:14 PM, Pierre Roudier <pierre.roudier@gmail.com>
wrote:

Bulent,

> Thank you for taking time and helping out. I attached a short version of
> data in (.xls) format. The file contains climate model results between 0
> and
> 800 BP, for every century, at a single weather station (Van-Turkey).

Thanks. This is much clearer :slight_smile:

> When you say "reformatting" do you actually mean organizing data for
> importing as
> a vector point?

Yes. Consider the following R session. I just exported your XLS data
sheet into CSV in LibreOffice.

> df <- read.csv('Van_Precipitation.csv')
> df
UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May
1 502500 4398240 1661 0 35.35559 34.10211 41.08517 54.96244 53.24204
2 502500 4398240 1661 -100 35.35434 34.10950 41.11511 54.98155 51.01671
3 502500 4398240 1661 -200 34.95282 34.74009 43.04097 56.12360 45.02433
4 502500 4398240 1661 -300 35.35213 34.12072 41.16004 55.01019 48.52734
5 502500 4398240 1661 -400 35.20788 34.33214 41.78447 55.40257 18.66510
6 502500 4398240 1661 -500 35.10954 34.49606 42.30229 55.70860 49.07919
7 502500 4398240 1661 -600 35.06167 34.58229 42.58510 55.86880 48.44339
8 502500 4398240 1661 -700 35.00383 34.69651 42.97523 56.08147 45.40806
9 502500 4398240 1661 -800 35.01142 34.68574 42.94451 56.06451 45.70229
Jun Jul Aug Sep Oct Nov Dec
Annual
1 17.729643 22.593137 0.7567056 11.717915 54.03994 41.26782 31.98477
398.8373
2 18.018884 22.797248 2.0428841 11.654079 53.91107 41.37436 31.99786
398.3736
3 10.662323 9.111732 4.5243246 4.254332 44.57887 48.04743 33.02368
368.0845
4 18.508462 23.089403 10.5560756 11.613746 53.75689 41.53623 32.01759
405.2488
5 21.111523 21.222403 19.7168880 10.175261 50.41907 43.68269 32.32056
384.0406
6 19.838364 3.158531 0.5505311 10.543923 47.83534 45.58733 32.58478
376.7945
7 15.610073 17.096418 8.2630408 12.148564 46.61419 46.60779 32.73860
395.6199
8 9.915775 8.735210 5.0867794 1.947931 45.02865 47.92435 32.96591
365.7697
9 10.045168 9.128862 5.3451146 2.768823 45.19434 47.83308 32.94591
367.6698

Then I'll use the reshape 2 to change the structure of df:

> library(reshape2)
> res <- melt(df, id.vars=c(1:4), variable.name='month',
> value.name='model_output')
> head(res, 15)
UTM_E UTM_N ELEV DATE month model_output
1 502500 4398240 1661 0 Jan 35.35559
2 502500 4398240 1661 -100 Jan 35.35434
3 502500 4398240 1661 -200 Jan 34.95282
4 502500 4398240 1661 -300 Jan 35.35213
5 502500 4398240 1661 -400 Jan 35.20788
6 502500 4398240 1661 -500 Jan 35.10954
7 502500 4398240 1661 -600 Jan 35.06167
8 502500 4398240 1661 -700 Jan 35.00383
9 502500 4398240 1661 -800 Jan 35.01142
10 502500 4398240 1661 0 Feb 34.10211
11 502500 4398240 1661 -100 Feb 34.10950
12 502500 4398240 1661 -200 Feb 34.74009
13 502500 4398240 1661 -300 Feb 34.12072
14 502500 4398240 1661 -400 Feb 34.33214
15 502500 4398240 1661 -500 Feb 34.49606

Here, I got the value of the result model for each month in each year,
for each station.

The following step is interpolation. This is not clear what you want
to do, but I guess what you want is to generate one raster layer per
month and per year.

I propose to use the plyr library, which allows you to split your data
in function of one or more variables (for us, month and year), and
apply a function to it.

Disclaimer: the following script is untested (that's just the general
idea) and should be launched from a R session WITHIN GRASS. Also, I
did not include the projection of the variables "grid" and "data".

# Generate an interpolation grid (the same for each raster to be
generated)
resolution <- 5 # your resolution
grid <- expand.grid(x = seq(min(df$UTM_E), max(df$UTM_E), by =
resolution), y = seq(min(df$UTM_N), max(df$UTM_N), by = resolution))
grid <- points2grid(grid)

# This is a function that generate a map from each "slice" of data
make_map <- function(data, grid){
# dependencies
require(sp)
require(rgdal)
require(spgrass6)
require(gstat) # for kriging
require(automap) # to automate kriging

# converting the data into a SpatialPointsDataFrame
coordinates(data) <- ~UTM_E + UTM_N

# INSERT YOUR FAVORITE INTERPOLATION METHOD HERE
res <- autoKrige(model_output ~ 1, data, grid)

# write it as a raster here
layer_name <- paste("model_output_", unique(data$month), "-",
unique(data$DATE), sep ="") # a name for the current layer
writeRAST6(res, layer_name, zcol="var1.pred", ignore.stderr=TRUE,
useGDAL=TRUE)

# return the result
return(res)
}

# Here is the plyr call - split the data, apply the function
"make_map", which interpolates the data, creates a layer name and
writes it in GRASS. Finally, it writes back everything in "rasters",
which is a list.
library(plyr)
rasters <- dlply(.data = res, .variables = .(month, DATE), .fun =
make_map, grid = grid)

Again, untested and probably buggy. Probably others can do better and
quicker with a bit of awk magic.

Hope this helps,

Pierre

> Cheers,
> Bulent
>
> On Sun, May 29, 2011 at 8:19 PM, Pierre Roudier
> <pierre.roudier@gmail.com>
> wrote:
>>
>> No you probably don't need to do so. What R allows you to do, is to
>> reformat your data frame so you do not need 400 files!
>>
>> Could you post a (short) overview of the structure of your data? From
>> what I understand, I'd use R and the plyr package to reformat your
>> data before attempting an import.
>>
>> Cheers,
>>
>> Pierre
>>
>> 2011/5/30 Bulent Arikan <bulent.arikan@gmail.com>:
>> > Hi Pierre,
>> > Thank you for answering. I use R too, so your reply was definitely
>> > helpful.
>> > I am a bit confused –and, I do not mean that I am asking for an
>> > answer
>> > from
>> > you when writing this– about the steps involving GRASS. So far I
>> > tried
>> > "r.in.ascii" and "r.in.bin" without success.Simply, the CSV file I
>> > have
>> > has
>> > 400 rows (each row represents one year) and I am not sure if I need
>> > to
>> > make
>> > individual CSV files for each year (means 400 files!!!) before
>> > importing
>> > into GRASS.
>> > Thank you for writing.
>> > Bulent
>> > On Sun, May 29, 2011 at 5:40 PM, Pierre Roudier
>> > <pierre.roudier@gmail.com>
>> > wrote:
>> >>
>> >> Hi Bulent,
>> >>
>> >> My first idea - but I may be biased as I come from the R world -
>> >> would
>> >> be to make the imports steps of your data from R. R has a XLS
>> >> reader(but I usually find it safer to first export your spreadsheet
>> >> to
>> >> a CSV file).
>> >>
>> >> Hope this helps,
>> >>
>> >> Pierre
>> >>
>> >> 2011/5/30 Bulent Arikan <bulent.arikan@gmail.com>:
>> >> > Dear List,
>> >> > I have an Excel file that contains precipitation data in 13
>> >> > columns
>> >> > (mean
>> >> > values for each month and annual average) from a region between
>> >> > certain
>> >> > years (in rows): eventually showing temporal changes in
>> >> > precipitation
>> >> > values. The data are taken at a weather station whose coordinates
>> >> > are
>> >> > known.
>> >> > I want to use this file in order to make an fcell raster map in
>> >> > GRASS.
>> >> > After
>> >> > some research I found few leads and mostly on interpolation of
>> >> > data.
>> >> > However, I am confused about the actual process turning the data
>> >> > in
>> >> > Excel
>> >> > into a raster map. I will greatly appreciate any suggestions as to
>> >> > which
>> >> > module to use and how to arrange the data.
>> >> >
>> >> > Thank you,
>> >> > --
>> >> > BÜLENT
>> >> >
>> >> > _______________________________________________
>> >> > grass-user mailing list
>> >> > grass-user@lists.osgeo.org
>> >> > http://lists.osgeo.org/mailman/listinfo/grass-user
>> >> >
>> >> >
>> >>
>> >>
>> >>
>> >> --
>> >> Scientist
>> >> Landcare Research, New Zealand
>> >
>> >
>> >
>> > --
>> > BÜLENT ARIKAN, PhD
>> > Postdoctoral Scholar
>> > Center for Social Dynamics and Complexity &
>> > School of Human Evolution and Social Change
>> > Arizona State University
>> > Tempe - AZ
>> > 85287-2402
>> >
>>
>>
>>
>> --
>> Scientist
>> Landcare Research, New Zealand
>
>
>
> --
> BÜLENT ARIKAN, PhD
> Postdoctoral Scholar
> Center for Social Dynamics and Complexity &
> School of Human Evolution and Social Change
> Arizona State University
> Tempe - AZ
> 85287-2402
>

--
Scientist
Landcare Research, New Zealand

--
BÜLENT ARIKAN, PhD
Postdoctoral Scholar
Center for Social Dynamics and Complexity &
School of Human Evolution and Social Change
Arizona State University
Tempe - AZ
85287-2402

--
Scientist
Landcare Research, New Zealand

6

Bulent wrote:

  UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May
1 502500 4398240 1661 0 35.35559 34.10211 41.08517 54.96244 53.24204
>> 2 502500 4398240 1661 -100 35.35434 34.10950 41.11511 54.98155 51.01671
3 502500 4398240 1661 -200 34.95282 34.74009 43.04097 56.12360 45.02433
4 502500 4398240 1661 -300 35.35213 34.12072 41.16004 55.01019 48.52734
[...]

to import that directly into GRASS, save as .csv out of Excel,
then:

v.in.ascii in=datafile.csv out=temperature_years \
   fs=tab x=2 y=3 skip=1 \
   columns='cat_id integer, utm_e double precision, utm_n double precision, elev integer, date_bp integer, temp_f_jan double precision, temp_f_feb double precision, [...], temp_f_dec double precision'

(fill in mar-nov, I didn't feel like typing it in..)

note "date_bp" is used instead of "date", to avoid the SQL
reserved-word conflict. for a quick try leave off the column
naming option and it will auto-detect the column types.

note set fs= space,tab,comma, etc as you exported from the
spreadsheet. Do not save using fixed width columns.
AFAIK v.in.ascii doesn't respect "quoting" around text strings.

Empty, NULL, nan, or "*" values for floating point no-data may be
problematic; full support for importing those is a work in
progress. For now if you have any try leaving them empty.

v.surf.rst does a very nice job as long as your starting points
are not highly clumpy (then you'll get segmentation artifacts and
have to start tweaking the interp coefficients). For the sample
data above it should be ok.

Pierre, you may be interested in looking at the v.kridge(.py)
module, behind the scenes it calls R from GRASS to do the
interpolation. Requires python-rpy2 & the GRASS<->R cran package
to be installed, and a number of other support packages it will
tell you about. (I'm having a little trouble starting it today..
but the method is sound)

Hamish

7

Bulent wrote:
> UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May
> 1 502500 4398240 1661 0 35.35559 34.10211 41.08517 54.96244 53.24204
> >> 2 502500 4398240 1661 -100 35.35434 34.10950 41.11511 54.98155 51.01671
> 3 502500 4398240 1661 -200 34.95282 34.74009 43.04097 56.12360 45.02433
> 4 502500 4398240 1661 -300 35.35213 34.12072 41.16004 55.01019 48.52734
> [...]

Hamish:

to import that directly into GRASS, save as .csv out of
Excel, then:

v.in.ascii in=datafile.csv out=temperature_years \

...

ah, sorry I didn't notice the e,n were the same for each year.
better to split it first using the tools you are familiar with
(thus R).

nevermind..,
Hamish

ps - If points are on a regular grid, r.in.xyz is an
alternative to v.in.ascii.

8

Hi Hamish,

Pierre, you may be interested in looking at the v.kridge(.py)
module, behind the scenes it calls R from GRASS to do the
interpolation. Requires python-rpy2 & the GRASS<->R cran package
to be installed, and a number of other support packages it will
tell you about. (I'm having a little trouble starting it today..
but the method is sound)

I am aware of that module (it was a GSoC project last year if I
remember well), even though I haven't used it in GRASS (yet). I guess
the limit is actually the fact that R would load the entire dataset
into memory - this can be a big headache with LiDAR data.

ah, sorry I didn't notice the e,n were the same for each year.
better to split it first using the tools you are familiar with
(thus R).

There actually more than one station in Bullent's data set, but the
difficulty that one raster layer must be interpolated for each month
of each year that has been modelled. Thus a bit of data tweaking in R.
No idea if awk/sed or any of those magic unix tools I don't know much
about could do the trick and avoid the call to R?

Pierre

Hamish

--
Scientist
Landcare Research, New Zealand

9

Pierre wrote:

There actually more than one station in Bullent's data set,
but the difficulty that one raster layer must be interpolated for
each month of each year that has been modelled. Thus a bit of data
tweaking in R.
No idea if awk/sed or any of those magic unix tools I don't
know much about could do the trick and avoid the call to R?

There's very little those power tools can't do to text files with a few
lines of code, but if one is already an expert in R, Python, Matlab, or
some other swiss army knife-of-a-scripting language, for one-offs or
prototypes you might as well use the tool you're most familiar with.
Many roads to get to Rome.

But perhaps, in this case & if the dataset is < ~3 million points, GRASS's
tools are the cleanest of all. As the original fit into a spreadsheet
I'd expect GRASS would have no problem with it.
Consider:

#make some test data, all at the same point in space
GRASS> cat << EOF > testpts
#cat,e,n,year_bp,temp_F
1,10,10,-100,32.1
2,10,10,-200,33.2
3,10,10,-300,34.3
EOF

#import them
GRASS> v.in.ascii in=testpts out=testpts cat=1 x=2 y=3 fs=, \
  columns='cat integer, easting double precision, northing double precision, year_bp integer, temp_f double precision'

#show the database
GRASS> v.db.select testpts
cat|easting|northing|year_bp|temp_f
1|10|10|-100|32.1
2|10|10|-200|33.2
3|10|10|-300|34.3

#select only the points matching a certain year into a new map
GRASS> v.extract in=testpts out=testpts_200bp where='year_bp = -200'

#show the database for the extraction
GRASS> v.db.select testpts_200bp
cat|easting|northing|year_bp|temp_f
2|10|10|-200|33.2

#interpolate from that extraction
GRASS> v.surf.rst in=testpts_200bp zcolumn=temp_f # ...

or maybe even skip the v.extract step and do the SQL query from the
original import vector map directly in v.surf.rst:

GRASS> v.surf.rst in=testpts_200bp zcolumn=temp_f where='year_bp = -200' # ...
?

of course you'll need a database for that so it won't work well for
millions and millions of data points (ie LIDAR)- awk or some custom C
program (etc) would then be needed to do the split efficiently. probably
python would do ok too.

I note that v.to.db has a where= option; I don't know about Excel but
OpenOffice can save directly to a DBF file. But again, many roads..

Hamish