hello,
because doing it by hand or with already existing GRASS commands ist too
tedious, Elshad Shirinov programmed a small module for me which fullfills
this purpose.
The problem was, that coordinates which has been taken from a corner or center
of a sampling plot represent in fact a 50x50m or 50x200m plot/transect.
This module (r.buffer.rect) creates an area around the coordinate with defined
x and y side length and with defined position of the input coordinate.
e.g.: you sampled a 50x50m plot but you took the coordinate (x) of the
bottom-right corner, hence you need the actual area of your plot (1), that is
what r.buffer.rect does:
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 x
hope this module will help other people as well, Martin
P.S.: I will add it to the GRASS add-on Wiki in the near future, when we
setup the *.html manual page.
Usage:
r.buffer.rect input=string output=string keyval=value x=value y=value
alignment=string [title="phrase"] [--overwrite] [--verbose] [--quiet]
Flags:
--o Allow output files to overwrite existing files
--v Verbose module output
--q Quiet module output
Parameters:
input Name of existing raster file
output Name for the output raster file
keyval Value of relevant pixels in the input raster
x Width of the area
y Height of the area
alignment Alignment of the pixel relative to the rectangle
options: center,top-left,top-right,bottom-left,bottom-right
title Title of the output raster file
On Montag, 16. Juni 2008 16:18:49 Martin Wegmann wrote:
> Hello,
>
> I am looking for a way to create a square buffer for a point vector so that
> each point is in the end the centroid of a square with defined side length.
>
> I though about r.mapcalc statements to grow north-south and west-east and
> then somehow create an area from this square. But I have no idea how to do
> it without a lot of hand-work.
>
> Any idea how to achieve that? Or hints where to find relevant documents.
>
> Martin
>
> _______________________________________________
> grass-user mailing list
> [email protected]
> http://lists.osgeo.org/mailman/listinfo/grass-user
#include "local_proto.h"
/*
Creation of rectangular area around a point raster
by Elshad Shirinov.
*/
struct alignment
{
char *name; /* method name */
char *text; /* menu display - full description */
int index;
};
static struct alignment alignments[] =
{
{"center", "Key pixel will be the center of the buffer.", 0},
{"top-left", "Key pixel will be top-left of the buffer.", 1},
{"top-right", "Key pixel will be the top-right of the buffer.", 2},
{"bottom-left", "Key pixel will be the bottom-left of the buffer.", 3},
{"bottom-right", "Key pixel will be the bottom-right of the buffer.", 4},
{0,0,0}
};
int main (int argc, char *argv[])
{
/* input */
char *oldname, *oldmapset;
/* input */
char *newname, *newmapset;
/* in and out file pointers */
int in_fd;
int out_fd;
/* parameters */
int keyval;
int x, y;
int align;
int sx, sy;
/* maps */
CELL *map, *newmap;
/* other parameters */
char *title;
/* helper variables */
int row, col;
CELL *result;
char* str;
int n, i;
RASTER_MAP_TYPE map_type;
struct Cell_head ch, window;
struct GModule *module;
struct
{
struct Option *input, *output;
struct Option *keyval, *x, *y;
struct Option *alignment;
struct Option *title;
} parm;
G_gisinit (argv[0]);
module = G_define_module();
module->keywords = _("raster");
module->description =
_("Creates rectangular area around input point raster.");
parm.input = G_define_option() ;
parm.input->key = "input" ;
parm.input->type = TYPE_STRING ;
parm.input->required = YES ;
parm.input->gisprompt = "old,cell,raster" ;
parm.input->description= _("Name of existing raster file") ;
parm.output = G_define_option() ;
parm.output->key = "output" ;
parm.output->type = TYPE_STRING ;
parm.output->required = YES ;
parm.output->gisprompt = "new,cell,raster,output" ;
parm.output->description= _("Name for the output raster file") ;
parm.keyval = G_define_option() ;
parm.keyval->key = "keyval" ;
parm.keyval->type = TYPE_INTEGER ;
parm.keyval->required = YES ;
parm.keyval->description= _("Value of relevant pixels in the input raster") ;
parm.x = G_define_option() ;
parm.x->key = "x" ;
parm.x->type = TYPE_INTEGER ;
parm.x->required = YES ;
parm.x->description= _("Width of the area") ;
parm.y = G_define_option() ;
parm.y->key = "y" ;
parm.y->type = TYPE_INTEGER ;
parm.y->required = YES ;
parm.y->description= _("Height of the area") ;
parm.alignment = G_define_option() ;
parm.alignment->key = "alignment" ;
parm.alignment->type = TYPE_STRING ;
parm.alignment->required = YES ;
str = parm.alignment->options = G_malloc(1024);
for (n = 0; alignments[n].name; n++)
{
if (n)
strcat (str, ",");
else
*str = 0;
strcat (str, alignments[n].name);
}
parm.alignment->description= _("Alignment of the key pixel relative to the new rectangular area") ;
parm.title = G_define_option() ;
parm.title->key = "title" ;
parm.title->key_desc = "\"phrase\"" ;
parm.title->type = TYPE_STRING ;
parm.title->required = NO ;
parm.title->description= _("Title of the output raster file") ;
if (G_parser(argc,argv))
exit(EXIT_FAILURE);
/* get name of input file */
oldname = parm.input->answer;
/* test input files existence */
if( (oldmapset = G_find_cell2(oldname,"")) == NULL )
{
G_warning ( _("%s: <%s> raster file not found\n"), G_program_name(), oldname);
G_usage();
exit(EXIT_FAILURE);
}
/* check if the new file name is correct */
newname = parm.output->answer;
if (G_legal_filename(newname) < 0)
{
G_warning("%s: <%s> illegal file name\n", G_program_name(), newname);
exit(EXIT_FAILURE);
}
newmapset = G_mapset();
/* read keyval */
sscanf(parm.keyval->answer, "%d", &keyval);
/* read x */
sscanf(parm.x->answer, "%d", &x);
/* read y */
sscanf(parm.y->answer, "%d", &y);
/* align with input */
if (G_get_cellhd(oldname, oldmapset, &ch) < 0) {
G_fatal_error ( _("Couldn't read CELL_HEAD structure from file <%s> in mapset <%s>\n"),
oldname, oldmapset);
exit(EXIT_FAILURE);
}
G_get_window(&window);
G_align_window(&window, &ch);
G_set_window(&window);
/* get size */
sx = G_window_cols();
sy = G_window_rows();
/* find alignment */
for (n = 0; (str = alignments[n].name); n++)
if (strcmp(str, parm.alignment->answer) == 0)
break;
if (str) {
align = alignments[n].index;
} else {
G_warning (_("<%s=%s> unknown %s"),
parm.alignment->key, parm.alignment->answer, parm.alignment->key);
G_usage();
exit(EXIT_FAILURE);
}
/* allocate map buffers */
map = (CELL*) G_malloc (sx * sy * sizeof(CELL));
newmap = (CELL*) G_malloc(sx * sy * sizeof(CELL));
result = G_allocate_c_raster_buf();
/* fill newmap with null */
G_set_c_null_value(newmap, sx * sy);
/* open map */
if ((in_fd = G_open_cell_old (oldname, oldmapset)) < 0)
{
G_fatal_error ( _("can't open cell file <%s> in mapset %s\n"), oldname, oldmapset);
G_usage();
exit(EXIT_FAILURE);
}
/* read map */
G_message("Reading map file:\n");
for (row = 0; row < sy; row++) {
G_get_c_raster_row (in_fd, map + row * sx, row);
G_percent (row + 1, sy, 1);
}
/* create buffers */
for(row = 0; row < sy; row++) {
for(col = 0; col < sx; col++) {
if(map[row * sx + col] == keyval) {
set_buffer(newmap, col, row, x, y, sx, sy, align);
}
}
}
/* close map */
G_close_cell(in_fd);
/* write the output file */
G_message("Writing output ... ");
/* open cell file */
if ((in_fd = G_open_cell_old (oldname, oldmapset)) < 0)
{
}
/* open new cell file */
out_fd = G_open_raster_new (newname, CELL_TYPE);
if (out_fd < 0) {
G_fatal_error("Can't create new cell file <%s> in mapset %s", newname, newmapset);
exit(EXIT_FAILURE);
}
/* write output */
for(row = 0; row < sy; row++) {
G_set_c_null_value(result, sx);
for(col = 0; col < sx; col++) {
result[col] = newmap[row * sx + col];
}
G_put_c_raster_row(out_fd, result);
G_percent (row + 1, sy, 1);
}
/* close new file */
G_close_cell(out_fd);
/* free allocated resources */
G_free(map);
G_free(newmap);
G_free(result);
exit(EXIT_SUCCESS);
}
MODULE_TOPDIR = ../..
PGM = r.buffer.rect
LIBES = $(STATSLIB) $(GISLIB)
DEPENDENCIES = $(STATSDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
default: cmd
#include "local_proto.h"
void set_buffer(CELL* buffer, int x, int y, int width, int height, int sx, int sy, int align) {
int l, r, t, b;
int i, j;
int dx, dy;
switch(align) {
case 0: /* center */
dx = width / 2;
dy = height / 2;
l = x - dx;
r = x + dx;
t = y - dy;
b = y + dy;
break;
case 1: /* top-left */
l = x;
r = x + width - 1;
t = y;
b = y + height - 1;
break;
case 2: /* top-right */
l = x - width + 1;
r = x;
t = y;
b = y + height - 1;
break;
case 3: /* bottom-left */
l = x;
r = x + width - 1;
t = y - height + 1;
b = y;
break;
case 4: /* bottom-right */
l = x - width + 1;
r = x;
t = y - height + 1;
b = y;
break;
default:
l = t = 0;
r = b = 1;
break;
}
l = l < 0 ? 0 : l;
r = r >= sx ? sx - 1 : r;
t = t < 0 ? 0 : t;
b = b >= sy ? sy - 1 : b;
/* fill buffer */
for(j = t; j <= b; j++) {
for(i = l; i <= r; i++) {
buffer[j * sx + i] = 1;
}
}
}
#ifndef LOCAL_PROTO_H
#define LOCAL_PROTO_H
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/stats.h>
#include <math.h>
#include <time.h>
/* buffer.c */
void set_buffer(CELL* buffer, int x, int y, int width, int height, int sx, int sy, int align);
#endif /* LOCAL_PROTO_H */
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user
