@Filename: focalcurve.c @Author: Logicon Geodynamics; Carol Gebhardt, based on an algorithm developed by: Jacek Blaszczynski, Physical Scientist, BLM NARSC and described in: Blaszczynski, J., 1997. Landform Characterization with GIS, PE&RS February 1997. @Non Portable Feature: None @Assumptions: None @Inputs: Ingrid File Name outgrid File Name Average slope grid File Name Negative average slope grid File Name Positive average slope grid File Name Negative count grid File Name Positive count grid File Name Kernel file Neighborhood shape Neighborhood cell width Neighborhood cell height Neighborhood outer radius Neighborhood inner radius Neighborhood quadrent direction @Outputs: Focalcurve grid Average slope grid Negative average slope grid Positive average slope grid Negative count grid Positive count grid @Description: Each cell of the input grid is processed using the specified neighborhood. For each neighborhood cell calculate the rise (value of center cell minus value of neighbor cell) and run (horizontal distance between the center cell and the neighbor cell). The slope is rise over run. Create the output grids from the accumulated slope values. ----------------------------------------------------------------- ------------ */ /* ----------------------------------------------------------------- ------------ INCLUDE FILES: ----------------------------------------------------------------- --------- */ #include #include #include #include #include #include #include /* ----------------------------------------------------------------- ------------ EXTERNAL DECLARATIONS: ----------------------------------------------------------------- --------- */ /* ----------------------------------------------------------------- ------------ DEFINES AND CONSTANTS: ----------------------------------------------------------------- --------- */ #define STATUS_OK 0 #define STATUS_ERROR 1 FILE *in_file; /* Input grid */ FILE *out_file; /* Output focal curve grid */ FILE *ave_file; /* Output average slope grid */ FILE *neg_ave_file; /* Output negative average ratio grid */ FILE *pos_ave_file; /* Output positive average ratio grid */ FILE *neg_count_file; /* Output negative count ratio grid */ FILE *pos_count_file; /* Output positive count ratio grid */ FILE *kernel_file; /* Input kernel file */ /* ----------------------------------------------------------------- ------------ STRUCTURES: ----------------------------------------------------------------- --------- */ /* * Command Line Parameters structure */ typedef struct _sParms { char ingrid [256]; /* Ingrid File Name */ char focalcurve_grid [256]; /* outgrid File Name */ char ave_slope_grid [256]; /* average slope grid File Name */ char neg_ave_grid [256]; /* negative average gird File Name */ char pos_ave_grid [256]; /* positive average gird File Name */ char neg_count_grid [256]; /* negative count gird File Name */ char pos_count_grid [256]; /* positive count gird File Name */ char kernel_file [256]; /* kernel file */ int shape; /* neighborhood shape */ long width; /* neighborhood cell width */ long height; /* neighborhood cell height */ float outer_radius; /* neighborhood outer radius */ float inner_radius; /* neighborhood inner radius */ int direction; /* neighborhood quadrent direction */ int min_x; /* minimum x neighborhood offset */ int max_x; /* maximum x neighborhood offset */ int min_y; /* minimum y neighborhood offset */ int max_y; /* maximum y neighborhood offset */ int ncols; /* number of columns in input grid */ int nrows; /* number of rows in input grid */ float xllcorner; /* x of lower left corner */ float yllcorner; /* y of lower left corner */ float cellsize; /* grid cell size */ float NODATA_value; /* value for no data */ } Parms; /* ** Neighborhood structure */ typedef struct _sNeigh { int row_num; /* neighborhood row number */ float *row_value; /* pointer to neighborhood row */ } Neigh; /* ** Kernel structure */ typedef struct _sKernel { float *kernel_value; /* kernel value */ } Kernel; /* ----------------------------------------------------------------- ------------ INTERNAL DECLARATIONS: ----------------------------------------------------------------- --------- */ /* ----------------------------------------------------------------- ------------ PRIVATE FUNCTION PROTOTYPES: ----------------------------------------------------------------- --------- */ /* ================================================================ OpenFiles */ /* ----------------------------------------------------------------- ------------ @Function: OpenFiles @Description: Open the input and output files and create output headers @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long OpenFiles ( /*RETURN - Completion Status */ Parms * parms); /*IN/OUT - input parameters */ /* ================================================================ GetKernel */ /* ----------------------------------------------------------------- ------------ @Function: GetKernel @Description: Get the kernel file and read into an array @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long GetKernel ( /*RETURN - Completion Status */ Parms * parms, /*IN - input parameters */ Kernel * kernel) /*OUT - kernel array */ ; /* ============================================================== ProcessGrid */ /* ----------------------------------------------------------------- ------------ @Function: ProcessGrid @Description: Process each of the grid cells to create new grids @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long ProcessGrid ( /*RETURN - Completion Status */ Parms * parms, /*IN - input parameters */ Kernel * kernel) /*IN - kernel array */ ; /* ========================================================== SetNeighborhood */ /* ----------------------------------------------------------------- ------------ @Function: SetNeighborhood @Description: Set the neighborhood boundaries @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long SetNeighborhood ( /*RETURN - Completion Status */ Parms * parms) /*IN/OUT - input parameters */ ; /* ====================================================== GetCommandLineParms */ /* ----------------------------------------------------------------- ------------ @Function: GetCommandLineParms @Description: Get and Validate the Command Line Parameters @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long GetCommandLineParms ( int argc, /*IN -Number of Arguments on command line */ char ** argv, /*IN -Command line arguments */ Parms * parms) /* OUT-Ptr to parameter structure */ ; /* ======================================================== focalcurve-main()*/ /* ----------------------------------------------------------------- ------------ @Function: focalcurve-main() @Description: @Notes: @Parameters: see below ----------------------------------------------------------------- ------------ */ int main ( /* RETURN-Completion Status: 0=Success; 1=Error */ int argc, /* N/A -Number of command line arguments (computed by system, not input by user */ char ** argv ) /*IN -Command line arguments ( i.e. usage ) 1 2 */ { Parms s_parms; /* Parameter structure */ Parms *parms=&s_parms; /* Ptr to Parameter struct */ Kernel s_kernel; /* Kernel structure */ Kernel *kernel = &s_kernel; /* Ptr to kernel array */ long status; /* Status variable */ /* ** Verify Command Line */ status = GetCommandLineParms( argc, argv, parms ); if ( status != STATUS_OK ) { return( status ); } /* ** Open the input and output files */ status = OpenFiles( parms ); if ( status != STATUS_OK ) { return( status ); } /* ** Set neighborhood boundaries */ status = SetNeighborhood( parms ); if ( status != STATUS_OK ) { return( status ); } /* ** Get kernel file */ status = GetKernel( parms, kernel ); if ( status != STATUS_OK ) { return( status ); } /* ** Process the input cells to create output grids */ status = ProcessGrid( parms, kernel ); if ( status != STATUS_OK ) { return( status ); } return( STATUS_OK ); } /* ----------------------------------------------------------------- ------------ PRIVATE FUNCTION BODIES: ----------------------------------------------------------------- --------- */ /* ============================================================== ProcessGrid */ /* ----------------------------------------------------------------- ------------ @Function: ProcessGrid @Description: Process each of the grid cells to create new grids @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long ProcessGrid ( Parms * parms, /*IN - input parameters */ Kernel * kernel) /*IN - kernel array */ { long status; /* Status variable */ long status_percent; /* Percent of rows completed */ long old_status; /* Previous percent of rows completed */ long current_col; /* Current grid column */ long current_row; /* Current grid row */ long neigh_col; /* Current neighborhood column */ long neigh_row; /* Current neighborhood row */ long current_row_index; /* Array index of current grid row */ long neigh_row_index; /* Array index of neighborhood grid row */ float center_cell; /* Current center cell value */ float neigh_cell; /* Neighborhood cell value */ long offset_x; /* Distance of neighbor from center */ long offset_y; /* Distance of neighbor from center */ Neigh grid_row[100]; /* Array of neighborhood rows */ long row_index; /* Index for neighborhood rows */ long col_index; /* Index for neighborhood columns */ long num_neighbor_cells; /* Number of cells in neighborhood */ long kernel_index; /* Kernel location that matches current */ /* cells */ float total_neg_slope; /* Total negative slope values */ float total_pos_slope; /* Total positive slope values */ long num_neg; /* Number of negative cells */ long num_pos; /* Number of positive cells */ float total_slope; /* Total slope values */ float total_abs_slope; /* Total absolute slope values */ float rise; /* Rise value between two cells */ float focalcurve_cell; /* Output cell value */ float neg_ave_cell; /* Output cell value */ float pos_ave_cell; /* Output cell value */ float neg_count_cell; /* Output cell value */ float pos_count_cell; /* Output cell value */ float ave_slope_cell; /* Output cell value */ float run; /* Distance between cells */ float slope; /* Rise over run */ /* ** Allocate space for neighborhood rows */ for (row_index = 0; row_index < parms->height; row_index++) { grid_row[row_index].row_value = (float *) malloc (sizeof(float) * parms->ncols); if(grid_row[row_index].row_value == 0) { return ( STATUS_ERROR); } } /* ** Initialize the neighborhood */ for (row_index = 0; row_index < parms->height; row_index++) { grid_row[row_index].row_num = row_index; for (col_index = 0; col_index < parms->ncols; col_index++) { fscanf(in_file,"%f ", &grid_row[row_index].row_value[col_index]); } fscanf(in_file,"\n"); } /* ** Start status meter */ status_percent = 0; printf("Focalcurvature: Cells Processed %%%d\n", status_percent); /* ** Loop through each input cell and caluclate the output cell ** from the neighborhood. */ for (current_row = 0; current_row < parms->nrows; current_row++) { /* ** Find the current row index (the neighborhood array index ** where the current row can be found) */ for (row_index = 0; row_index < parms->height; row_index++) { if (grid_row[row_index].row_num == current_row) { current_row_index = row_index; } } /* ** Loop through the cells in the current row */ for (current_col = 0; current_col < parms->ncols; current_col++) { /* ** Initialize the variables */ total_neg_slope = 0; total_pos_slope = 0; num_neg = 0; num_pos = 0; total_slope = 0; total_abs_slope = 0; num_neighbor_cells = 0; focalcurve_cell = 0; ave_slope_cell = 0; neg_count_cell = 0; pos_count_cell = 0; neg_ave_cell = 0; pos_ave_cell = 0; kernel_index = 0; /* ** Get the value of the center cell */ center_cell = grid_row[current_row_index].row_value[current_col]; if(center_cell != parms->NODATA_value) { /* ** Loop through the neighborhood cells */ for (neigh_row = current_row + parms->min_y; neigh_row <= current_row + parms->max_y; neigh_row++) { /* ** Make sure neighborhood cell is on grid */ if (neigh_row > -1 && neigh_row < parms->nrows) { /* ** Find the neighborhood row index (the neighborhood array index ** where the neighborhood row can be found) */ for (row_index = 0; row_index < parms->height; row_index++) { if (grid_row[row_index].row_num == neigh_row) { neigh_row_index = row_index; } } for (neigh_col = current_col + parms->min_x; neigh_col <= current_col + parms->max_x; neigh_col++) { /* ** Make sure neighborhood cell is on grid */ if (neigh_col > -1 && neigh_col < parms->ncols) { neigh_cell = grid_row[neigh_row_index].row_value[neigh_col]; /* ** If center cell then skip it, ** or if kernel = NODATA skip it ** or if neighbor = NODATA skip it */ if ((neigh_row == current_row && neigh_col == current_col) || kernel->kernel_value[kernel_index] == parms->NODATA_value || neigh_cell == parms->NODATA_value) { kernel_index++; } else { neigh_cell = neigh_cell * kernel->kernel_value[kernel_index]; rise = center_cell - neigh_cell; offset_x = fabs(current_col - neigh_col); offset_y = fabs(current_row - neigh_row); run = sqrt((parms->cellsize * offset_y) * (parms->cellsize * offset_y) + (parms->cellsize * offset_x) * (parms->cellsize * offset_x)); num_neighbor_cells ++; if(run != 0) { slope = (rise/run); total_slope = total_slope + slope; total_abs_slope = total_abs_slope + fabs(slope); if (slope < 0) { num_neg ++; total_neg_slope = total_neg_slope + slope; } else { num_pos++; total_pos_slope = total_pos_slope + slope; } } else { printf("Divide by Zero error\n"); } kernel_index++; } } else { /* ** Neighborhood cell is not on grid. ** Increment kernel */ kernel_index++; } } /* end loop through neighborhood columns */ } /* end if row in neighborhood */ else { /* ** Increment by whole row if row off grid */ kernel_index = kernel_index + parms->width; } } /* end loop through the neighborhood rows */ /* ** Create new grid cell values */ if (num_neighbor_cells > 0) { focalcurve_cell = total_slope / num_neighbor_cells; fprintf(out_file, "%f ", focalcurve_cell); } else { fprintf(out_file, "%f ", 0); } if(strncmp(parms->ave_slope_grid, "NO_FILE", 7)) { if(num_neighbor_cells != 0) { ave_slope_cell = total_abs_slope / num_neighbor_cells; fprintf(ave_file, "%f ", ave_slope_cell); } else { fprintf(ave_file, "%f ", 0); } } if(strncmp(parms->pos_ave_grid, "NO_FILE", 7)) { if(num_pos == 0) { pos_ave_cell = 0; } else { pos_ave_cell = total_pos_slope / num_pos; } fprintf(pos_ave_file, "%f ", pos_ave_cell); } if(strncmp(parms->neg_ave_grid, "NO_FILE", 7)) { if(num_neg == 0) { neg_ave_cell = 0; } else { neg_ave_cell = total_neg_slope / num_neg; } fprintf(neg_ave_file, "%f ", neg_ave_cell); } if(strncmp(parms->neg_count_grid, "NO_FILE", 7)) { fprintf(neg_count_file, "%d ", num_neg); } if(strncmp(parms->pos_count_grid, "NO_FILE", 7)) { fprintf(pos_count_file, "%d ", num_pos); } } else { /* ** Center is NODATA so output is NODATA */ fprintf(out_file, "%f ", (float)parms->NODATA_value); if(strncmp(parms->ave_slope_grid, "NO_FILE", 7)) { fprintf(ave_file, "%f ", (float)parms->NODATA_value); } if(strncmp(parms->neg_count_grid, "NO_FILE", 7)) { fprintf(neg_count_file, "%f ", (float)parms->NODATA_value); } if(strncmp(parms->pos_count_grid, "NO_FILE", 7)) { fprintf(pos_count_file, "%f ", (float)parms->NODATA_value); } if(strncmp(parms->neg_ave_grid, "NO_FILE", 7)) { fprintf(neg_ave_file, "%f ", (float)parms->NODATA_value); } if(strncmp(parms->pos_ave_grid, "NO_FILE", 7)) { fprintf(pos_ave_file, "%f ", (float)parms->NODATA_value); } } } /* end loop through each column */ fprintf(out_file, "\n"); if(strncmp(parms->ave_slope_grid, "NO_FILE", 7)) { fprintf(ave_file, "\n"); } if(strncmp(parms->neg_count_grid, "NO_FILE", 7)) { fprintf(neg_count_file, "\n"); } if(strncmp(parms->pos_count_grid, "NO_FILE", 7)) { fprintf(pos_count_file, "\n"); } if(strncmp(parms->neg_ave_grid, "NO_FILE", 7)) { fprintf(neg_ave_file, "\n"); } if(strncmp(parms->pos_ave_grid, "NO_FILE", 7)) { fprintf(pos_ave_file, "\n"); } /* ** Move neighborhood down one row (if not end of file) */ if (current_row + parms->max_y + 1 < parms->nrows) { for (row_index = 0; row_index < parms->height; row_index++) { if (grid_row[row_index].row_num == current_row + parms->min_y) { grid_row[row_index].row_num = current_row + parms->max_y + 1; for (col_index = 0; col_index < parms->ncols; col_index++) { fscanf(in_file,"%f ", &grid_row[row_index].row_value[col_index]); } fscanf(in_file,"\n"); } } } /* if more rows */ /* ** Print status of processing */ old_status = status_percent; status_percent = (int)((((float)current_row + 1) / (float)parms->nrows) * 100); if(status_percent != old_status && status_percent % 25 == 0) { printf("Focalcurvature: Cells Processed %%%d\n", status_percent); } } /* loop through rows */ /* ** Close files */ fclose(in_file); fclose(out_file); if(strncmp(parms->ave_slope_grid, "NO_FILE", 7)) { fclose(ave_file); } if(strncmp(parms->neg_count_grid, "NO_FILE", 7)) { fclose(neg_count_file); } if(strncmp(parms->pos_count_grid, "NO_FILE", 7)) { fclose(pos_count_file); } if(strncmp(parms->neg_ave_grid, "NO_FILE", 7)) { fclose(neg_ave_file); } if(strncmp(parms->pos_ave_grid, "NO_FILE", 7)) { fclose(pos_ave_file); } /* ** Free space for (row_index = 0; row_index < parms->height; row_index++) { free (grid_row[row_index].row_value); } */ free (kernel->kernel_value); /* ** Everything completed OK */ return( STATUS_OK ); } /* ====================================================== GetCommandLineParms */ /* ----------------------------------------------------------------- ------------ @Function: GetCommandLineParms @Description: Get the arguments from the command line @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long GetCommandLineParms ( /*RETURN-Completion Status */ int argc, /*IN - Number of Arguments in Command Line */ char ** argv, /*IN - Command Line Arguments */ Parms * parms) /*OUT - Ptr to parameter struct */ { long status; /* Status returned from functions */ /* ** Make sure there are enough parameters */ if ( argc != 15 ) { fprintf(stderr, "Insufficient parameters for focalcurve\n"); return( STATUS_ERROR ); } /* ** Process the Command Line */ strcpy(parms->ingrid, argv[1]); strcpy(parms->focalcurve_grid, argv[2]); strcpy(parms->ave_slope_grid, argv[3]); strcpy(parms->neg_ave_grid, argv[4]); strcpy(parms->pos_ave_grid, argv[5]); strcpy(parms->neg_count_grid, argv[6]); strcpy(parms->pos_count_grid, argv[7]); parms->shape = atoi(argv[8]); parms->width = atoi(argv[9]); parms->height = atoi(argv[10]); parms->outer_radius = atof(argv[11]); parms->inner_radius = atof(argv[12]); parms->direction = atoi(argv[13]); strcpy(parms->kernel_file, argv[14]); return( STATUS_OK ); } /* ================================================================ GetKernel */ /* ----------------------------------------------------------------- ------------ @Function: GetKernel @Description: The kernel array has two uses: 1) If irregular or weighted, the user will supply a kernel file which is read into the kernel array. 2) If no kernel file is supplied, the kernel array is used to mask the shape desired by the user. NODATA is put into the array for cells that are not to be concidered part of the neighborhood. @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long GetKernel ( /*RETURN-Completion Status */ Parms * parms, /*IN - input parameters */ Kernel * kernel) /*OUT - kernel array */ { long status; /* Status returned from functions */ int num_kernels; /* Number of kernels */ int kernel_index; /* Kernel index */ int current_y; /* Current kernel y location */ int current_x; /* Current kernel x location */ float run; /* Distance between cells */ Kernel s_kernel; /* Kernel structure */ Kernel * my_kernel = &s_kernel;/* Ptr to kernel structure */ if (strncmp(parms->kernel_file, "NO_FILE", 7)) { /* ** Open the kernel File. */ if( (kernel_file = fopen( parms->kernel_file, "r" )) == NULL) { printf("Cannot open '%s' \n", parms->kernel_file ); return( STATUS_ERROR ); } /* ** Read the header */ fscanf(kernel_file, "%d %d\n", &parms->width, &parms->height); num_kernels = parms->width * parms->height; /* ** Create the space for the kernel array */ my_kernel->kernel_value = (float *) malloc (sizeof(float)*num_kernels); if (my_kernel->kernel_value == 0) { return ( STATUS_ERROR); } /* ** If shape is WEIGHT then read in kernel as written */ if (parms->shape == 6) { kernel_index = 0; /* ** Read file into array */ for (current_y = 0; current_y < parms->height; current_y++) { for (current_x = 0; current_x < parms->width; current_x++) { fscanf(kernel_file, "%f ", &my_kernel->kernel_value[kernel_index]); kernel_index++; } fscanf(kernel_file,"\n"); } } else { kernel_index = 0; /* ** If IRREGULAR and kernel value is zero, change to the ** NODATA_value. If kernel value is non-zero, value is 1. */ for (current_y = 0; current_y < parms->height; current_y++) { for (current_x = 0; current_x < parms->width; current_x++) { fscanf(kernel_file, "%f ", &my_kernel->kernel_value[kernel_index]); if(my_kernel->kernel_value[kernel_index] == 0) { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; } else { my_kernel->kernel_value[kernel_index] = 1; } kernel_index++; } fscanf(kernel_file,"\n"); } } /* ** Close file */ fclose(kernel_file); } else /* no kernel file */ { num_kernels = parms->width * parms->height; /* ** Create the space for the kernel array */ if (NULL==(my_kernel->kernel_value = (float *)malloc(num_kernels * sizeof(float)))) { return ( STATUS_ERROR); } kernel_index = 0; /* ** If rectangle, all kernel values are 1 */ if (parms->shape == 1) { for (kernel_index = 0; kernel_index < parms->width * parms->height; kernel_index++) { my_kernel->kernel_value[kernel_index] = 1.0; } } /* ** If circle, kernel values are 1 if within outer radius */ if (parms->shape == 2) { for (current_y = parms->min_y; current_y <= parms->max_y; current_y++) { for (current_x = parms->min_x; current_x <= parms->max_x; current_x++) { run = sqrt(current_x * current_x + current_y * current_y); if (run <= parms->outer_radius) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } } } /* ** If annulus, kernel values are 1 if smaller than outer radius ** and bigger than inner radius */ if (parms->shape == 3) { for (current_y = parms->min_y; current_y <= parms->max_y; current_y++) { for (current_x = parms->min_x; current_x <= parms->max_x; current_x++) { run = sqrt(current_x * current_x + current_y * current_y); if (run <= parms->outer_radius && run >= parms->inner_radius) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } } } /* ** If wedge, kernel values are 1 if smaller than outer radius ** and within the specified quarter circle */ if (parms->shape == 4) { for (current_y = parms->min_y; current_y <= parms->max_y; current_y++) { for (current_x = parms->min_x; current_x <= parms->max_x; current_x++) { run = sqrt(current_x * current_x + current_y * current_y); if (run <= parms->outer_radius) { if(parms->direction == 1) /* east */ { if(current_x > 0) { if(fabs(atan(current_y / current_x)) <= atan(1)) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else if(parms->direction == 2) { if(current_y >= 0 && current_x != 0) { if(fabs(atan(current_y / current_x)) >= atan(1) && fabs(atan(current_y / current_x)) <= atan(-1)) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else if(parms->direction == 3) { if(current_x < 0) { if(fabs(atan(current_y / current_x)) <= atan(1)) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else if(parms->direction == 4) { if(current_y <= 0 && current_x != 0) { if(fabs(atan(current_y / current_x)) >= fabs(atan(1)) && fabs(atan(current_y / current_x)) <= fabs(atan(-1))) { my_kernel->kernel_value[kernel_index] = 1.0; kernel_index++; } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } } else { my_kernel->kernel_value[kernel_index] = parms->NODATA_value; kernel_index++; } } } } } kernel->kernel_value = my_kernel->kernel_value; /* KERNEL DEBUG PRINT for (kernel_index = 0; kernel_index < num_kernels; kernel_index++) { printf("kernel %f \n", kernel->kernel_value[kernel_index]); } */ return( STATUS_OK ); } /* ================================================================ OpenFiles */ /* ----------------------------------------------------------------- ------------ @Function: OpenFiles @Description: Open the input and output files and create output headers @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long OpenFiles ( /*RETURN-Completion Status */ Parms * parms) /*IN - input parameters */ { long status; /* Status returned from functions */ long index; /* String index */ /* ** Open the ingrid File. */ if( (in_file = fopen( parms->ingrid, "r" )) == NULL) { /* ** If can't open file, try converting to lower case first. */ for(index = 0; index <= strlen(parms->ingrid); index++) { parms->ingrid[index] = tolower(parms->ingrid[index]); } if( (in_file = fopen( parms->ingrid, "r" )) == NULL) { printf("Cannot open '%s' \n", parms->ingrid ); return( STATUS_ERROR ); } } /* ** Open the outgrid File. */ if( (out_file = fopen( parms->focalcurve_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->focalcurve_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fscanf(in_file,"ncols %d\n", &parms->ncols); fprintf(out_file,"ncols %d\n",parms->ncols); fscanf(in_file,"nrows %d\n", &parms->nrows); fprintf(out_file,"nrows %d\n",parms->nrows); fscanf(in_file,"xllcorner %f\n", &parms->xllcorner); fprintf(out_file,"xllcorner %.3f\n",parms->xllcorner); fscanf(in_file,"yllcorner %f\n", &parms->yllcorner); fprintf(out_file,"yllcorner %.3f\n",parms->yllcorner); fscanf(in_file,"cellsize %f\n", &parms->cellsize); fprintf(out_file,"cellsize %.3f\n",parms->cellsize); fscanf(in_file,"NODATA_value %f\n", &parms->NODATA_value); fprintf(out_file,"NODATA_value %.3f\n",parms->NODATA_value); /* ** Open the negative ratio File. */ if (strncmp(parms->neg_ave_grid, "NO_FILE", 7)) { if( (neg_ave_file = fopen( parms->neg_ave_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->neg_ave_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fprintf(neg_ave_file,"ncols %d\n",parms->ncols); fprintf(neg_ave_file,"nrows %d\n",parms->nrows); fprintf(neg_ave_file,"xllcorner %.3f\n",parms->xllcorner); fprintf(neg_ave_file,"yllcorner %.3f\n",parms->yllcorner); fprintf(neg_ave_file,"cellsize %.3f\n",parms->cellsize); fprintf(neg_ave_file,"NODATA_value %.3f\n",parms->NODATA_value); } /* ** Open the positive ratio File. */ if (strncmp(parms->pos_ave_grid, "NO_FILE", 7)) { if( (pos_ave_file = fopen( parms->pos_ave_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->pos_ave_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fprintf(pos_ave_file,"ncols %d\n",parms->ncols); fprintf(pos_ave_file,"nrows %d\n",parms->nrows); fprintf(pos_ave_file,"xllcorner %.3f\n",parms->xllcorner); fprintf(pos_ave_file,"yllcorner %.3f\n",parms->yllcorner); fprintf(pos_ave_file,"cellsize %.3f\n",parms->cellsize); fprintf(pos_ave_file,"NODATA_value %.3f\n",parms->NODATA_value); } /* ** Open the negative count File. */ if (strncmp(parms->neg_count_grid, "NO_FILE", 7)) { if( (neg_count_file = fopen( parms->neg_count_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->neg_count_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fprintf(neg_count_file,"ncols %d\n",parms->ncols); fprintf(neg_count_file,"nrows %d\n",parms->nrows); fprintf(neg_count_file,"xllcorner %.3f\n",parms->xllcorner); fprintf(neg_count_file,"yllcorner %.3f\n",parms->yllcorner); fprintf(neg_count_file,"cellsize %.3f\n",parms->cellsize); fprintf(neg_count_file,"NODATA_value %.3f\n",parms->NODATA_value); } /* ** Open the positive count File. */ if (strncmp(parms->pos_count_grid, "NO_FILE", 7)) { if( (pos_count_file = fopen( parms->pos_count_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->pos_count_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fprintf(pos_count_file,"ncols %d\n",parms->ncols); fprintf(pos_count_file,"nrows %d\n",parms->nrows); fprintf(pos_count_file,"xllcorner %.3f\n",parms->xllcorner); fprintf(pos_count_file,"yllcorner %.3f\n",parms->yllcorner); fprintf(pos_count_file,"cellsize %.3f\n",parms->cellsize); fprintf(pos_count_file,"NODATA_value %.3f\n",parms->NODATA_value); } /* ** Open the average slope File. */ if (strncmp(parms->ave_slope_grid, "NO_FILE", 7)) { if( (ave_file = fopen( parms->ave_slope_grid, "w" )) == NULL) { printf("Cannot open '%s' \n", parms->ave_slope_grid ); return( STATUS_ERROR ); } /* ** Create the output file headers */ fprintf(ave_file,"ncols %d\n",parms->ncols); fprintf(ave_file,"nrows %d\n",parms->nrows); fprintf(ave_file,"xllcorner %.3f\n",parms->xllcorner); fprintf(ave_file,"yllcorner %.3f\n",parms->yllcorner); fprintf(ave_file,"cellsize %.3f\n",parms->cellsize); fprintf(ave_file,"NODATA_value %.3f\n",parms->NODATA_value); } return( STATUS_OK ); } /* ========================================================== SetNeighborhood */ /* ----------------------------------------------------------------- ------------ @Function: SetNeighborhood @Description: Set the neighborhood boundaries @Notes: @Parameters: See below ----------------------------------------------------------------- ------------ */ static long SetNeighborhood ( /*RETURN - Completion Status */ Parms * parms) /*IN/OUT - input parameters */ { long status; /* Status returned from functions */ long x_offset; /* X offset from center */ long y_offset; /* Y offset from center */ /* ** Set the neighborhood boundaries */ x_offset = (int)((parms->width -1) / 2); y_offset = (int)((parms->height -1) / 2); parms->min_x = x_offset * (-1); parms->min_y = y_offset * (-1); if((parms->width -1) % 2 > 0) { parms->max_x = x_offset + 1; } else { parms->max_x = x_offset; } if((parms->height -1) % 2 > 0) { parms->max_y = y_offset + 1; } else { parms->max_y = y_offset; } return( STATUS_OK ); }