[GRASS-user] making fcell raster maps using climate data
Pierre Roudier
pierre.roudier at gmail.com
Mon May 30 02:14:38 EDT 2011
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 :)
> 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 at 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 at 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 at 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 at 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 at 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
More information about the grass-user
mailing list