[GRASSLIST:1546] Re: counting pixels

[Glynn Clements]

Uttam Kumar wrote:

>                   I have a raster map. It contains 720 rows x 960 columns 
> =691200 pixels.
> There are 32676 categories. I want to generate the category for each pixel; 
> somewhat in this manner. For one pixel the lat-long bound comes to be a grid 
> of 2 min,15 sec x 2 min, 15 sec. I want the output somewhat like this:
> -------------------------------------------------------------------
> pixel          North lat  South lat   East long   West long     category
> ------         --------   ---------   ----------  ----------    --------
> 
> 1,1           07:02:15N   07:00:00N   68:02:15E    68:00:00E       34
> 1,2           07:02:15N   07:00:00N   68:04:30E    68:02:15E       21
> ..            ........    ........     ........    ........        ..
> ..            ........    ........     ........    ........        ..
> 720,960       34:50:15     34:48:00    93:02:00    93:00:00        97
> 
> I am trying to first take the first row and traversing all the 960 columns, 
> then taking the second row and traversing through all the 960 column and 
> continuing the process till all the 720 rows and 960 columns are covered. I 
> think we can get a similar output by using g.region and ultimately setting 
> the region of each pixel by its lat long value (perhaps by a loop). Yet I 
> want all the data in the single file unlike generating a separate file for 
> each category. Since creating separate file for each pixel would result in 
> 691200 files....
> very difficult to read. Is it possible in GRASS ? Please suggest me some 
> solution.

Try the attached C program. Compile with e.g.:

gcc grass_dump.c -o grass_dump -I/usr/local/grass5/include -L/usr/local/grass5/lib -lgis

-- 
Glynn Clements <glynn.clements at virgin.net>

-------------- next part --------------

#include <stdio.h>
#include <stdlib.h>

#include "gis.h"

static void lat_format(double lat, char *buf)
{
	int d, m;
	double s;
	char h;

	G_lat_parts(lat, &d, &m, &s, &h);

	sprintf(buf, "%02d:%02d:%02d%c", d, m, (int)(s + 0.5), h);
}

static void lon_format(double lon, char *buf)
{
	int d, m;
	double s;
	char h;

	G_lon_parts(lon, &d, &m, &s, &h);

	sprintf(buf, "%03d:%02d:%02d%c", d, m, (int)(s + 0.5), h);
}

int main(int argc, char **argv)
{
	struct Cell_head window;
	char *name, *mapset;
	int fd;
	CELL *inbuf;
	int x, y;

	struct GModule *module;
	struct Option *input;

	G_gisinit(argv[0]);

	module = G_define_module();
	module->description =
		"Dump a GRASS raster file as a table";

	input = G_define_option();
	input->key		= "input";
	input->type		= TYPE_STRING;
	input->required		= YES;
	input->multiple		= NO;
	input->gisprompt	= "old,cell,raster";
	input->description	= "Raster file to be read.";

	if (G_parser(argc, argv))
		exit(1);

	name = input->answer;

	mapset = G_find_cell2(name, "");
	if (!mapset)
		G_fatal_error("Couldn't find raster file %s", name);

	fd = G_open_cell_old(name, mapset);
	if (fd < 0) 
		G_fatal_error("Couldn't open raster file %s@%s", name, mapset);

	inbuf = G_allocate_cell_buf();

	G_get_window(&window);

	printf("%-12s%-12s%-12s%-12s%-12s %s\n",
	       "pixel", "North lat", "South lat", "East long", "West long", "category");

	printf("%-12s%-12s%-12s%-12s%-12s %s\n",
	       "-----", "---------", "---------", "---------", "---------", "--------");

	for (y = 0; y < window.rows; y++)
	{
		double n = window.north - window.ns_res * y;
		double s = window.north - window.ns_res * (y + 1);
		char nbuf[12], sbuf[12];
		char yxbuf[12];

		lat_format(n, nbuf);
		lat_format(s, sbuf);

		G_get_c_raster_row(fd, inbuf, y);

		for (x = 0; x < window.cols; x++)
		{
			double w = window.west + window.ew_res * x;
			double e = window.west + window.ew_res * (x + 1);
			CELL cat = inbuf[x];
			char wbuf[12], ebuf[12];

			if (G_is_c_null_value(&cat))
				cat = -1;

			sprintf(yxbuf, "%d,%d", y, x);

			lon_format(w, wbuf);
			lon_format(e, ebuf);

			printf("%-12s%-12s%-12s%-12s%-12s %d\n",
				yxbuf, nbuf, sbuf, ebuf, wbuf, cat);
		}
	}

	G_close_cell(fd);

	G_free(inbuf);

	return 0;
}