'bugfix' for s.surf.idw
Job Spijker
spijker at geo.uu.nl
Thu Oct 28 10:51:09 EDT 1999
On Thu, 28 Oct 1999, Stephan Eickschen wrote:
> I think that many people will be interested in your code - so why don't
> you send it to the list - or directly to me - so it will be in the
> archive very soon? Then everybody wil be able to make a download... Did
> you here something from Markus or Baylor (well, the latter are not awake
> already, I guess :) )?
The code is attached to this mail, I don't know if everyone on the list
appreciates it to receive large chunks of code so that's why I didn't
sent it to the list (I'm on the list for a about two weeks so I'm not
familiar with the habbits of the list). As I said it's a quick and dirty
hack and for example it doesn't take care of masks nor is it validated but
I think if you don't use masks no errors will occur.
I didn't hear from Markus or Baylor, you where the first to contact me.
Grtz, Job
----------------------------------------------------------------------
Job Spijker
Faculty of Earth Sciences/ Department of Geochemistry
Utrecht University, The Netherlands
Budapestlaan 4 3584 CD Utrecht
T: 030-253 3264 F: 030-253 5030
----------------------------------------------------------------------
-------------- next part --------------
/* file changed 19991028 by Job Spijker (spijker at geo.uu.nl)
* changes are marked with 'hack' */
#include "gis.h"
int search_points = 12;
int npoints = 0;
int npoints_alloc = 0;
int nsearch;
struct Point
{
double north, east;
double z;
double dist;
};
struct Point *points = NULL;
struct Point *list;
int main(int argc, char *argv[])
{
int fd, maskfd;
CELL *cell, *mask;
DCELL *dcell; /*hack*/
struct Cell_head window;
int row, col;
double north, east;
double dx,dy;
double maxdist,dist;
double sum1, sum2;
double dummy; /*hack*/
int i,n,max;
void read_sites();
struct
{
struct Option *input, *npoints, *output;
} parm;
parm.input = G_define_option() ;
parm.input->key = "input" ;
parm.input->type = TYPE_STRING ;
parm.input->required = YES ;
parm.input->description= "Name of input sites map" ;
parm.input->gisprompt = "old,site_lists,sites" ;
parm.output = G_define_option() ;
parm.output->key = "output" ;
parm.output->type = TYPE_STRING ;
parm.output->required = YES;
parm.output->description= "Name of output raster map" ;
parm.output->gisprompt = "any,cell,raster" ;
parm.npoints = G_define_option() ;
parm.npoints->key = "npoints" ;
parm.npoints->key_desc = "count" ;
parm.npoints->type = TYPE_INTEGER ;
parm.npoints->required = NO ;
parm.npoints->description="Number of interpolation points";
parm.npoints->answer = "12";
G_gisinit(argv[0]);
if (G_parser(argc, argv))
exit(1);
if (G_legal_filename(parm.output->answer) < 0)
{
fprintf (stderr, "%s=%s - illegal name\n", parm.output->key, parm.output->answer);
exit(1);
}
if(sscanf(parm.npoints->answer,"%d", &search_points) != 1 || search_points<1)
{
fprintf (stderr, "%s=%s - illegal number of interpolation points\n",
parm.npoints->key, parm.npoints->answer);
G_usage();
exit(1);
}
list = (struct Point *) G_calloc (search_points, sizeof (struct Point));
/* read the elevation points from the input sites file */
read_sites (parm.input->answer);
if (npoints == 0)
{
fprintf (stderr, "%s: no data points found\n", G_program_name());
exit(1);
}
nsearch = npoints < search_points ? npoints : search_points;
/* get the window, allocate buffers, etc. */
G_get_set_window (&window);
dcell=G_allocate_d_raster_buf(); /*hack*/
cell = G_allocate_cell_buf();
if ((maskfd = G_maskfd()) >= 0)
mask = G_allocate_cell_buf();
else
mask = NULL;
/*fd = G_open_cell_new(parm.output->answer);*/
fd=G_open_raster_new(parm.output->answer,2); /*hack*/
if (fd < 0)
{
fprintf (stderr, "%s: can't create %s\n", G_program_name(), parm.output->answer);
exit(1);
}
fprintf (stderr, "Interpolating raster map <%s> ... %d rows ... ",
parm.output->answer, window.rows);
north = window.north + window.ns_res/2.0;
for (row = 0; row < window.rows; row++)
{
fprintf (stderr, "%-10d\b\b\b\b\b\b\b\b\b\b", window.rows-row);
if (mask)
{
if(G_get_map_row(maskfd, mask, row) < 0)
exit(1);
}
north -= window.ns_res;
east = window.west - window.ew_res/2.0;
for (col = 0; col < window.cols; col++)
{
east += window.ew_res;
/* don't interpolate outside of the mask */
if (mask && mask[col] == 0)
{
cell[col] = 0;
continue;
}
/* fill list with first nsearch points */
for (i = 0; i < nsearch ; i++)
{
dy = points[i].north - north;
dx = points[i].east - east;
list[i].dist = dy*dy + dx*dx;
list[i].z = points[i].z;
}
/* find the maximum distance */
maxdist = list[max=0].dist;
for (n = 1; n < nsearch; n++)
{
if (maxdist < list[n].dist)
maxdist = list[max=n].dist;
}
/* go thru rest of the points now */
for ( ; i < npoints; i++)
{
dy = points[i].north - north;
dx = points[i].east - east;
dist = dy*dy + dx*dx;
if (dist < maxdist)
{
/* replace the largest dist */
list[max].z = points[i].z;
list[max].dist = dist;
maxdist = list[max=0].dist;
for (n = 1; n < nsearch; n++)
{
if (maxdist < list[n].dist)
maxdist = list[max=n].dist;
}
}
}
/* interpolate */
sum1 = 0.0;
sum2 = 0.0;
for (n = 0; n < nsearch; n++)
{
if(dist = list[n].dist)
{
sum1 += list[n].z / dist;
sum2 += 1.0/dist;
}
else
{
sum1 = list[n].z;
sum2 = 1.0;
break;
}
}
dcell[col] = (DCELL) (sum1/sum2 + .5);
/* for debug purposes */
/*
dummy=((sum1/sum2)+0.5);
fprintf (stderr,"INTERPOLATE calculcated cel %i,%i : %f dummy :%f\n",row,col,dcell[col],dummy);
*/
}
/*G_put_map_row (fd, cell);*/
G_put_d_raster_row(fd,dcell); /*hack*/
}
G_close_cell(fd);
fprintf (stderr, "done \n");
exit(0);
}
void newpoint ( double z,double east,double north)
{
if (npoints_alloc <= npoints)
{
npoints_alloc += 128;
points = (struct Point *) G_realloc (points,
npoints_alloc * sizeof (struct Point));
}
points[npoints].north = north;
points[npoints].east = east;
points[npoints].z = z;
npoints++;
}
More information about the grass-user
mailing list