Re: [GRASS-user] making fcell raster maps using climate data
Hi Bulent, No problem, but please always copy the list, so that: - others can offer an alternative solution (as I said, I'm highly biased towards R) - our discussion can be searchable by everybody for further reference > I have spent a good portion of today working on this project and I got as > far as reorganizing the data. It seems that any step after this is beyond my > knowledge of GRASS, such as generating an interpolation grid (quoting you): > # Generate an interpolation grid (the same for each raster to be generated) > resolution <- 5 # your resolution This just sets the resolution of the grid I want to generate to 5m. > 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 quick and dirty way to generate a grid (raster) on which interpolate your points. If you got a polygon (shapefile) of your study area, you'd better use it to generate an interpolation grid. See ?spsample for that. > I understand that I should define the geographic extents –hence creating an > interpolation grid– of the data I will interpolate from this weather > station. However, I am not sure if the statements below are the commands > that I should type in the Terminal of Grass (they seem to be R commands): > [ 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) ] > and > [ make_map <- function(data, grid) ] This is R code, indeed. What you need to do is to open a R session from within GRASS. Just type R on your GRASS command line. Then you'll be in a R session that is in GRASS :) Make sure you check out the spgrass6 package to interface your R session with your GRASS session. HTH, Pierre -- Scientist Landcare Research, New Zealand ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
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 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
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 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
> Bulent wrote: > > UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May > > 1 502500 4398240 16610 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. ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
Bulent wrote: > UTM_E UTM_N ELEV DATE Jan Feb Mar Apr May > 1 502500 4398240 16610 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 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
> Pierre,
> I cannot thank you enough for all these detailed, step-by-step explanations.
No worries ;)
> 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
> 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 :)
>>
>> > 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
Re: [GRASS-user] making fcell raster maps using climate data
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 16610 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 JulAug 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 monthmodel_output
1 502500 4398240 16610 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 16610 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,
Pierr
Re: [GRASS-user] making fcell raster maps using climate data
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 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 > ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
You can convert the xls file to csv file. Next, add annotation to see it as a vector (points) layer in gdal/ogr (http://www.gdal.org/ogr/drv_csv.html). After importing the data to grass you can interpolate the data to obtain precipitation field. You can use built-in interpolation methods (for flat and homogeneous terrain) or you can use regression based on terrain characteristics (for annual and seasonal sums it works very well). Yet another option is use geostatistical methods. Pawel - "Bulent Arikan" napisał: > 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 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] making fcell raster maps using climate data
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 : > 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 ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[GRASS-user] making fcell raster maps using climate data
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
