[GRASS5] Re: [GRASSLIST:6443] d.profile: phenology curve plotting
Jachym Cepicky
jachym.cepicky at centrum.cz
Thu Apr 14 11:01:40 EDT 2005
hallo developers,
On Wed, Apr 13, 2005 at 01:58:30PM -0400, Martin Wegmann wrote:
> Hello,
>
> is there a function which plots the phenology of a time series?
>
> d.profile is doing it spatially but I would like to look at the temporal curve
> of various raster which are a time series.
>
> e.g. d.profile input=raster.jan, raster.feb, .... raster.dec and then select a
> pixel/area in the image and receive a time profile for each pixel which I
> selected.
>
> thanks in advance, Martin
>
> --
> Martin Wegmann
>
I thought, that r.profile would be the best for this problem.
Unfortunally it works only with one raster.
I tryed to write some patch, that makes r.profile to handel multiple raster
files.
before:
r.profile in=dem_int out=- profile=106.25,893.75,866.6,139.58 res=100
Using Resolution 100.000000
Output Format:
[Along Track Dist.(m)] [Elevation]
Approx. transect length 1070.936279 m.
0.000000 -249
100.000000 -260
200.000000 -222
300.000000 -164
400.000000 -133
500.000000 -176
600.000000 -78
700.000000 252
800.000000 230
900.000000 199
1000.000000 28
After:
r.profile in=dem,dem4_int,dem3 out=- profile=106.25,893.75,866.6,139.58 res=100
Using Resolution 100.000000
Output Format:
[Along Track Dist.(m)] [Elevation]
Approx. transect length 1070.936279 m.
0.000000 -249.258438 -344 -152.711473
100.000000 -260.661446 -280 -137.916396
200.000000 -222.767889 -43 125.295406
300.000000 -164.659004 294 93.851748
400.000000 -133.269469 600 112.628475
500.000000 -176.494114 732 132.492686
600.000000 -78.644163 685 129.169117
700.000000 252.142540 410 76.520597
800.000000 230.412563 18 7.827528
900.000000 199.797620 -116 -87.963766
1000.000000 28.419509 -304 -316.939340
I thing, it could be usefull change.
Patches are in attachment.
Do you thing, it could be useful change?
Jachym
--
Jachym Cepicky
e-mail: jachym.cepicky at centrum.cz
URL: http://les-ejk.cz
GPG: http://www.fle.czu.cz/~jachym/gnupg_public_key/
-------------- next part --------------
Index: main.c
===================================================================
RCS file: /home/grass/grassrepository/grass51/raster/r.profile/main.c,v
retrieving revision 2.1
diff -u -r2.1 main.c
--- main.c 15 Jan 2005 11:18:13 -0000 2.1
+++ main.c 14 Apr 2005 14:54:58 -0000
@@ -17,11 +17,15 @@
#include "raster.h"
#include "glocale.h"
+#ifndef MAXFILE
+#define MAXFILE 10
+#endif
+
static int move(int, int);
static int cont(int, int);
-int do_profile(double, double, double, double, char *, int, double, int, int, FILE *, char *);
+int do_profile(double, double, double, double, char *[], int, double, int fd[], int data_type[], FILE *, char *, int );
/* read_rast.c */
-int read_rast(double, double, double, int, int, RASTER_MAP_TYPE, FILE *, char *);
+int read_rast(double, double, double, int fd[], int, RASTER_MAP_TYPE data_type[], FILE *, char *, int);
static int min_range[5], max_range[5];
static int which_range;
@@ -30,15 +34,15 @@
int main(int argc, char *argv[])
{
- char *name, *outfile, *mapset, msg[256];
- int fd, projection;
+ char *name[MAXFILE], *outfile, *mapset, msg[256];
+ int fd[MAXFILE], projection;
FILE *fp;
int screen_x, screen_y, button;
double res;
char errbuf[256], *null_string;
int coords = 0, i, k = -1;
double e1, e2, n1, n2;
- RASTER_MAP_TYPE data_type;
+ RASTER_MAP_TYPE data_type[MAXFILE];
struct Cell_head window;
struct
{
@@ -48,6 +52,9 @@
parm;
struct GModule *module;
+ char **ptr;
+ int maxfile;
+
G_gisinit(argv[0]);
/* Set description */
@@ -59,7 +66,7 @@
parm.opt1->key = "input";
parm.opt1->type = TYPE_STRING;
parm.opt1->required = YES;
- parm.opt1->multiple = NO;
+ parm.opt1->multiple = YES;
parm.opt1->gisprompt = "old,cell,raster";
parm.opt1->description = _("Name of existing raster map");
@@ -111,6 +118,7 @@
G_fatal_error(msg);
}
+
null_string = parm.null_str->answer;
G_get_window(&window);
@@ -135,19 +143,34 @@
G_begin_distance_calculations();
/* Open Input File for reading */
+
+ ptr = parm.opt1->answers;
+
/* Get Input Name */
- name = parm.opt1->answer;
+ for (maxfile = 0; *ptr != NULL; maxfile++, ptr++) {
+ if (maxfile >= MAXFILE)
+ G_fatal_error("%s - too many files. only %d allowed",
+ G_program_name(), MAXFILE);
+
+ name[maxfile] = *ptr;
+ }
+
if (parm.g->answer)
coords = 1;
- /* Open Raster File */
- if (NULL == (mapset = G_find_cell2(name, ""))) {
- sprintf(msg, "Cannot find mapset for %s \n", name);
- G_fatal_error(msg);
+ /* Open Raster Files */
+ for (i = 0; i < maxfile; i++) {
+ if (NULL == (mapset = G_find_cell2(name[i], ""))) {
+ sprintf(msg, "Cannot find mapset for %s \n", name[i]);
+ G_fatal_error(msg);
+ }
}
- if (0 > (fd = G_open_cell_old(name, mapset))) {
- sprintf(msg, "Cannot open File %s\n", name);
- G_fatal_error(msg);
+
+ for (i = 0; i < maxfile; i++) {
+ if (0 > (fd[i] = G_open_cell_old(name[i], mapset))) {
+ sprintf(msg, "Cannot open File %s\n", name[i]);
+ G_fatal_error(msg);
+ }
}
/* Open ASCII file for output or stdout */
@@ -162,7 +185,8 @@
}
/* Get Raster Type */
- data_type = G_raster_map_type(name, mapset);
+ for (i = 0; i < maxfile; i++)
+ data_type[i] = G_raster_map_type(name[i], mapset);
/* Done with file */
/* Get coords */
@@ -217,7 +241,7 @@
G_plot_line(e1, n1, e2, n2);
/* Get profile info */
- do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type, fp, null_string);
+ do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type, fp, null_string, maxfile);
n1 = n2;
e1 = e2;
@@ -240,7 +264,7 @@
e2 = e1;
n2 = n1;
/* Get profile info */
- do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type, fp, null_string);
+ do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type, fp, null_string, maxfile);
}
else {
for (i = 0; i <= k - 2; i += 2) {
@@ -250,13 +274,17 @@
sscanf(parm.profile->answers[i + 3], "%lf", &n2);
/* Get profile info */
do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type,
- fp, null_string);
+ fp, null_string, maxfile);
}
}
}
- G_close_cell(fd);
+ putchar('\n');
+ for (i = 0; i < maxfile; i++) {
+ G_close_cell(fd[i]);
+ }
+
fclose(fp);
return 0;
@@ -264,8 +292,7 @@
/* Calculate the Profile Now */
/* Establish parameters */
-int do_profile(double e1, double e2, double n1, double n2, char *name, int coords,
- double res, int fd, int data_type, FILE * fp, char *null_string)
+int do_profile(double e1, double e2, double n1, double n2, char *name[], int coords, double res, int fd[], int data_type[], FILE * fp, char *null_string, int maxfile)
{
float rows, cols, LEN;
double Y, X, AZI;
@@ -274,14 +301,14 @@
rows = n1 - n2;
LEN = G_distance(e1, n1, e2, n2);
- fprintf(stderr, "Approx. transect length %f m.\n", LEN);
+ fprintf(stderr, "Approx. transect length %f m.", LEN);
/* Calculate Azimuth of Line */
if (rows == 0 && cols == 0) {
/* Special case for no movement */
e = e1;
n = n1;
- read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+ read_rast(e, n, dist, fd, coords, data_type, fp, null_string, maxfile);
}
if (rows >= 0 && cols < 0) {
@@ -297,7 +324,7 @@
dist -= G_distance(e, n, e1, n1);
}
for (e = e1, n = n1; e < e2 || n > n2; e += X, n -= Y) {
- read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+ read_rast(e, n, dist, fd, coords, data_type, fp, null_string, maxfile);
/* d+=res; */
dist += G_distance(e - X, n + Y, e, n);
}
@@ -319,7 +346,7 @@
*/
}
for (e = e1, n = n1; e < e2 || n < n2; e += X, n += Y) {
- read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+ read_rast(e, n, dist, fd, coords, data_type, fp, null_string, maxfile);
/* d+=res; */
dist += G_distance(e - X, n - Y, e, n);
}
@@ -338,7 +365,7 @@
dist -= G_distance(e, n, e1, n1);
}
for (e = e1, n = n1; e > e2 || n > n2; e -= X, n -= Y) {
- read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+ read_rast(e, n, dist, fd, coords, data_type, fp, null_string, maxfile);
/* d+=res; */
dist += G_distance(e + X, n + Y, e, n);
}
@@ -357,7 +384,7 @@
dist -= G_distance(e, n, e1, n1);
}
for (e = e1, n = n1; e > e2 || n < n2; e -= X, n += Y) {
- read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+ read_rast(e, n, dist, fd, coords, data_type, fp, null_string, maxfile);
/* d+=res; */
dist += G_distance(e + X, n - Y, e, n);
}
-------------- next part --------------
Index: read_rast.c
===================================================================
RCS file: /home/grass/grassrepository/grass51/raster/r.profile/read_rast.c,v
retrieving revision 2.0
diff -u -r2.0 read_rast.c
--- read_rast.c 9 Nov 2004 13:51:36 -0000 2.0
+++ read_rast.c 14 Apr 2005 14:55:04 -0000
@@ -17,8 +17,8 @@
#include "gis.h"
#include "display.h"
-int read_rast(double east, double north, double dist, int fd, int coords,
- RASTER_MAP_TYPE data_type, FILE * fp, char *null_string)
+int read_rast(double east, double north, double dist, int fd[], int coords,
+ RASTER_MAP_TYPE data_type[], FILE * fp, char *null_string, int maxfile)
{
int row, col, nrows, ncols;
struct Cell_head window;
@@ -26,6 +26,7 @@
char buf[1024] = "";
FCELL *fcell;
DCELL *dcell;
+ int i;
G_get_window(&window);
nrows = window.rows;
@@ -37,51 +38,77 @@
if (row < 0)
G_fatal_error("Coordinate request outsite current region settings");
- if (data_type == CELL_TYPE) {
- cell = G_allocate_c_raster_buf();
- if (G_get_c_raster_row(fd, cell, row) < 0)
- exit(1);
-
- if (G_is_c_null_value(&cell[col]))
- sprintf(buf, null_string);
- else
- sprintf(buf, "%d", cell[col]);
-
- if (coords == 1)
- fprintf(fp, "%f %f %f %s\n", east, north, dist, buf);
- else
- fprintf(fp, "%f %s\n", dist, buf);
- }
-
- if (data_type == FCELL_TYPE) {
- fcell = G_allocate_f_raster_buf();
- if (G_get_f_raster_row(fd, fcell, row) < 0)
- exit(1);
-
- if (G_is_f_null_value(&fcell[col]))
- sprintf(buf, null_string);
- else
- sprintf(buf, "%f", fcell[col]);
-
- if (coords == 1)
- fprintf(fp, "%f %f %f %s\n", east, north, dist, buf);
- else
- fprintf(fp, "%f %s\n", dist, buf);
- }
-
- if (data_type == DCELL_TYPE) {
- dcell = G_allocate_d_raster_buf();
- if (G_get_d_raster_row(fd, dcell, row) < 0)
- exit(1);
- if (G_is_d_null_value(&dcell[col]))
- sprintf(buf, null_string);
- else
- sprintf(buf, "%f", dcell[col]);
-
- if (coords == 1)
- fprintf(fp, "%f %f %f %s\n", east, north, dist, buf);
- else
- fprintf(fp, "%f %s\n", dist, buf);
+ /* for each raster file */
+ for (i = 0; i < maxfile; i++) {
+
+ if (data_type[i] == CELL_TYPE) {
+ cell = G_allocate_c_raster_buf();
+ if (G_get_c_raster_row(fd[i], cell, row) < 0)
+ exit(1);
+
+ if (G_is_c_null_value(&cell[col]))
+ sprintf(buf, null_string);
+ else
+ sprintf(buf, "%d ", cell[col]);
+
+ if (coords == 1) {
+ if (i == 0) /* print coordninates and
+ distance only for the first raster */
+ fprintf(fp, "\n%f %f %f ", east, north, dist);
+ fprintf(fp, "%s ", buf);
+ }
+ else {
+ if (i == 0) /* print coordinates and
+ distance only for the first raster */
+ fprintf(fp, "\n%f ", dist);
+ fprintf(fp, "%s", buf);
+ }
+ }
+
+ if (data_type[i] == FCELL_TYPE) {
+ fcell = G_allocate_f_raster_buf();
+ if (G_get_f_raster_row(fd[i], fcell, row) < 0)
+ exit(1);
+
+ if (G_is_f_null_value(&fcell[col]))
+ sprintf(buf, null_string);
+ else
+ sprintf(buf, "%f", fcell[col]);
+
+ if (coords == 1) {
+ if (i == 0)
+ fprintf(fp, "\n%f %f %f ", east, north, dist);
+ fprintf(fp, "%s ", buf);
+ }
+ else {
+ if (i == 0)
+ fprintf(fp, "\n%f", dist);
+ fprintf(fp, "%s", buf);
+ }
+ }
+
+ if (data_type[i] == DCELL_TYPE) {
+ dcell = G_allocate_d_raster_buf();
+
+ if (G_get_d_raster_row(fd[i], dcell, row) < 0)
+ exit(1);
+
+ if (G_is_d_null_value(&dcell[col]))
+ sprintf(buf, null_string);
+ else
+ sprintf(buf, "%f", dcell[col]);
+
+ if (coords == 1) {
+ if (i == 0)
+ fprintf(fp, "\n%f %f %f ", east, north, dist);
+ fprintf(fp, "%s\t", buf);
+ }
+ else {
+ if (i == 0)
+ fprintf(fp, "\n%f ", dist);
+ fprintf(fp, "%s ", buf);
+ }
+ }
}
return 0;
More information about the grass-dev
mailing list