001/* 002 * 003 * Copyright (C) 2001 Numeric Solutions. All Rights Reserved. 004 * 005 * Permission is hereby granted, without written agreement and without license 006 * or royalty fees, to use, copy, modify, and distribute this software and its 007 * documentation for any purpose, provided that the above copyright notice 008 * appears in all copies of this software. 009 * 010 * Contact information: Numeric Solutions support@numericsolutions.com, or: 011 * 012 * http://www.numericsolutions.com 013 * 014 */ 015 016package com.numericsolutions.geomodeltools; 017 018import java.io.BufferedReader; 019import java.io.FileOutputStream; 020import java.io.PrintStream; 021import java.io.StringReader; 022import java.util.StringTokenizer; 023import java.util.Vector; 024 025public class invdist_power_isosearch2d { 026 027 private int NX = 0; 028 private int NY = 0; 029 030 /** 031 * Method that takes as input a name of a file and writes the file contents 032 * to a vector and then makes the vector and number of vector elements 033 * access to the public 034 * 035 * @param fileName 036 * name of the file that whose contents should be written to a 037 * vector 038 * @return -- the contents of the file as a Vector 039 */ 040 private Vector StringVectorizer(String data) { 041 Vector localVector = new Vector(); 042 try { 043 BufferedReader in = new BufferedReader(new StringReader(data)); 044 String s; 045 while ((s = in.readLine()) != null) { 046 localVector.addElement(s); 047 } 048 } catch (Exception e) { 049 System.out.println("failed: " + e.getMessage()); 050 } 051 return (localVector); 052 } 053 054 /** 055 * Creates a Raster grid in memory with the following structure: 056 * 057 * x, y, p, i, j 058 * 059 * @param xmin 060 * float -- the minumum x value 061 * @param xmax 062 * float -- the maximum x value 063 * @param ymin 064 * float -- the minimum y value 065 * @param ymax 066 * float -- the maximum y value 067 * @param dx 068 * float -- the distance between cells in the x direction 069 * @param dy 070 * float -- the distance between cells in the y direction 071 * @param p 072 * float -- the default grid value - this will be replaced by the 073 * interpolation process 074 * @return Vector -- the grid vector, with the afore mentioned structure 075 */ 076 private Vector getTemplateGrid(float xmin, float xmax, float ymin, 077 float ymax, float dx, float dy, float p) { 078 079 int nx = (int) ((xmax - xmin) / dx); // gonna loose precision here 080 int ny = (int) ((ymax - ymin) / dy); 081 082 // set the instance vars 083 this.NX = nx; 084 this.NY = ny; 085 086 Vector out = new Vector(); 087 088 for (int i = 1; i < ny + 1; i++) { 089 float y = ymin + (dy * i); 090 for (int ii = 1; ii < nx + 1; ii++) { 091 float x = xmin + (dx * ii); 092 String row = x + " " + y + " " + p + " " + ii + " " + i; 093 out.addElement(row); 094 } 095 } 096 return (out); 097 } 098 099 /** 100 * This is the interpolation function for this class, it does an Inverse 101 * Distance Weighted (IDW) Interpolation, using a user-defined weight 102 * coefficient, and serach radius. 103 * 104 * @param gridVec 105 * Vector -- a raster grid with an x, y, p, i, j structure 106 * 107 * @param data 108 * String -- the datafile 109 * 110 * @param outFile 111 * String -- the output file 112 * 113 * @param weight 114 * float -- the weight to be used for the interpolation, the 115 * smaller this value (or closer to 1) the smoother the grid, the 116 * larger the value the better the fit to the data. Typically, 117 * for geological data a value of 2 is used and for biological 118 * data 1.2 is commonly used. 119 * 120 * @param searchRadius 121 * float -- the lateral search direction -- if a grid cell has no 122 * corresponding data within the search no interpolation on that 123 * cell is done 124 * 125 * @return Vector -- the raster grid as a vector with the structure: x, y, 126 * value, i, j 127 */ 128 public Vector getInverseDistance2dGrid(Vector gridVec, String data, 129 String outFile, float weight, float searchRadius) { 130 131 // System.err.println("weight --> " + weight + 132 // "\nsearch radius --> " + searchRadius); 133 134 Vector outVec = new Vector(); 135 136 try { 137 Vector dataVec = this.StringVectorizer(data); 138 PrintStream out = new PrintStream(new FileOutputStream(outFile)); 139 float power = (weight); 140 double searchradius = searchRadius; 141 double point_estimate; 142 int datalines = dataVec.size(); 143 int gridlines = gridVec.size(); 144 String datavals[] = new String[datalines]; 145 String gridvals[] = new String[gridlines]; 146 147 dataVec.toArray(datavals); 148 gridVec.toArray(gridvals); 149 150 double distvals[] = new double[datalines]; 151 152 int cnt = 0; 153 while (cnt < gridlines) { 154 StringTokenizer t = new StringTokenizer(gridvals[cnt], " \t"); 155 double cur_x = (new Double(t.nextToken())).doubleValue(); 156 double cur_y = (new Double(t.nextToken())).doubleValue(); 157 double cur_p = (new Double(t.nextToken())).doubleValue(); 158 int i = Integer.parseInt(t.nextToken()); 159 int j = Integer.parseInt(t.nextToken()); 160 161 int data_cnt = 0; 162 while (data_cnt < datalines) { 163 StringTokenizer tdata = new StringTokenizer( 164 datavals[data_cnt], " \t"); 165 double cur_datax = (new Double(tdata.nextToken())) 166 .doubleValue(); 167 double cur_datay = (new Double(tdata.nextToken())) 168 .doubleValue(); 169 double cur_datap = (new Double(tdata.nextToken())) 170 .doubleValue(); 171 172 double dy = cur_datay - cur_y; 173 double dx = cur_datax - cur_x; 174 if (dx == 0.0) { 175 dx = 0.000001; 176 } 177 /* 178 * if dx, the numerator is zero make it some very small 179 * number 180 */ 181 if (dy == 0.0) { 182 dy = 0.000001; 183 } 184 /* this is to keep dy from equal zero at 90 degrees */ 185 double dy_dx = dy / dx; 186 double angle_radians = Math.atan(dy_dx); 187 double angle_degrees = angle_radians / (Math.PI / 180); 188 /* 189 * convert zero angle to some small number - this should / 190 * not ever happen because of the above logic if statements 191 */ 192 193 double dist = dy / (Math.sin(angle_radians)); 194 /* load the array with the distnces and angles */ 195 distvals[data_cnt] = dist; 196 data_cnt++; 197 } 198 199 /* 200 * next block will calculate the weights to be used for the 201 * averaging 202 */ 203 int array_cnt = 0; 204 double invdist_sum = 0.0; 205 double weights[] = new double[datalines]; /* 206 * this is the weights 207 * array 208 */ 209 while (array_cnt < datalines) { 210 if (distvals[array_cnt] < 0) { 211 distvals[array_cnt] = distvals[array_cnt] * (-1); 212 } 213 if (distvals[array_cnt] <= searchradius) { 214 invdist_sum = invdist_sum 215 + (1 / Math.pow(distvals[array_cnt], power)); 216 } 217 array_cnt++; 218 } 219 220 /* 221 * next bit divides each inverse distance by the sum of the 222 * inverses 223 */ 224 array_cnt = 0; 225 while (array_cnt < datalines) { 226 /* 227 * if statemant checks to see that the values are less than 228 * the serach 229 */ 230 if (distvals[array_cnt] <= searchradius) { 231 weights[array_cnt] = (1 / Math.pow(distvals[array_cnt], 232 power)) 233 / invdist_sum; 234 } 235 if (distvals[array_cnt] > searchradius) { 236 weights[array_cnt] = 0; 237 } 238 array_cnt++; 239 } 240 241 /* 242 * next block multiplies the weights by the values to calculate 243 * the estimates 244 */ 245 array_cnt = 0; 246 point_estimate = 0; 247 while (array_cnt < datalines) { 248 StringTokenizer datapoints = new StringTokenizer( 249 datavals[array_cnt], " \t"); 250 String buf1 = datapoints.nextToken(); 251 buf1 = datapoints.nextToken(); 252 double cur_dataz = (new Double(datapoints.nextToken())) 253 .doubleValue(); 254 point_estimate = point_estimate 255 + (weights[array_cnt] * cur_dataz); 256 array_cnt++; 257 } 258 259 outVec.addElement(cur_x + " " + cur_y + " " + point_estimate 260 + " " + i + " " + j); 261 out.println(cur_x + " " + cur_y + " " + point_estimate + " " 262 + i + " " + j); 263 cnt++; 264 } 265 out.close(); 266 } catch (Exception e) { 267 e.printStackTrace(); 268 } 269 270 return outVec; 271 } 272 273 /** 274 * this is the main method that is used to create the arc grid file from the 275 * fileName 276 * 277 * @param xyzFile 278 * -- this is the file that contains the xyz values delimeted by 279 * a space 280 */ 281 private StringBuffer writeArcGrid(Vector gridVec, double xmin, double xmax, 282 double ymin, double ymax, int rows, int cols, double resolution, 283 double nullValue) { 284 285 StringBuffer sb = new StringBuffer(); 286 try { 287 288 double x = xmin; 289 double y = ymax; 290 Vector z_vec = new Vector(); 291 Vector row_vec = new Vector(); 292 Vector col_vec = new Vector(); 293 Vector matrix = new Vector(); 294 295 // poopulate the vectors that are to be used for searching 296 for (int k = 0; k < gridVec.size(); k++) { 297 298 z_vec.addElement(getColumnToken((String) gridVec.elementAt(k), 299 3)); 300 row_vec.addElement(getColumnToken( 301 (String) gridVec.elementAt(k), 4)); 302 col_vec.addElement(getColumnToken( 303 (String) gridVec.elementAt(k), 5)); 304 305 } 306 307 // the header 308 /** 309 * System.out.println("ncols " + rows); // notice that these 310 * are switched System.out.println("nrows " + cols); 311 * System.out.println("xllcorner " + xmin); 312 * System.out.println("yllcorner " + ymin); 313 * System.out.println("cellsize " + resolution); 314 * System.out.println("NODATA_value " + nullValue); 315 **/ 316 317 sb.append("ncols " + rows + "\n"); // notice that these are 318 // switched 319 sb.append("nrows " + cols + "\n"); 320 sb.append("xllcorner " + xmin + "\n"); 321 sb.append("yllcorner " + ymin + "\n"); 322 sb.append("cellsize " + resolution + "\n"); 323 sb.append("NODATA_value " + nullValue + "\n"); 324 325 String longLine = ""; 326 for (int i = 1; i <= cols; i++) { 327 for (int j = 1; j <= rows; j++) { 328 for (int k = 0; k < gridVec.size(); k++) { 329 // if this is the node that we're searching for 330 if (i == Integer.parseInt(((String) col_vec 331 .elementAt(k))) 332 && j == Integer.parseInt(((String) row_vec 333 .elementAt(k)))) { 334 String z = (String) z_vec.elementAt(k); 335 if (j == 1) { 336 longLine = z; 337 } else { 338 longLine = longLine + " " + z; 339 } 340 break; 341 } 342 } 343 } 344 // System.out.println(longLine); 345 matrix.addElement(longLine); 346 } 347 for (int i = 1; i <= matrix.size(); i++) { 348 sb.append(matrix.elementAt(matrix.size() - i) + "\n"); 349 // System.out.println(matrix.elementAt(matrix.size() - i)); 350 } 351 } catch (Exception e) { 352 e.printStackTrace(); 353 } 354 return sb; 355 } 356 357 /** 358 * utility method that returns a column value from a string that contains 359 * multiple tokens 360 * 361 * @param string 362 * -- the string to tokenize 363 * @param tokenPosition 364 * -- the location in the string to get 365 */ 366 public String getColumnToken(String s, int tokenPosition) { 367 String token = "null"; 368 if (s != null) { 369 // System.out.println( "not null"); 370 StringTokenizer t = new StringTokenizer(s.trim(), " \r\n\t"); 371 int i = 1; 372 while (i <= tokenPosition) { 373 if (t.hasMoreTokens()) { 374 token = t.nextToken(); 375 i++; 376 } else { 377 token = "null"; 378 i++; 379 } 380 } 381 return (token); 382 } else { 383 return (token); 384 } 385 } 386 387 /** 388 * returns the application usage as a String 389 * 390 * @return String -- the string describing how to use the application 391 */ 392 private static String getUsage() { 393 StringBuffer sb = new StringBuffer(); 394 sb.append("\n\nusage: $java invdist_power_isosearch2d params \n"); 395 sb.append(" where params: "); 396 sb.append(" xmin, xmax, ymin, ymax, dx, dy, null value, datafile, " 397 + "output file, weight coeficient, search radius\n\n"); 398 399 sb 400 .append("[xmin] --> the minimum x location describing the output grid\n"); 401 sb 402 .append("[xmax] --> the maximum x location describing the output grid\n"); 403 sb 404 .append("[ymin] --> the minimum y location describing the output grid\n"); 405 sb 406 .append("[ymax] --> the maximum y location describing the output grid\n"); 407 sb 408 .append("[dx] --> the spacing between the grid cells in the x direction\n"); 409 sb 410 .append("[dy] --> the spacing between the grid cells in the y direction\n"); 411 sb.append("[nullvalue] --> the default null value\n"); 412 413 sb 414 .append("[datafile] --> the data file used for the interpolation -- this \n"); 415 sb.append("file needs to have a structure of: x y property like: \n\n"); 416 sb.append("34.94365 -119.7953333 772.8 \n" 417 + "34.94365 -119.2928333 769.8 \n" 418 + "34.9438333 -119.9943333 918.1 \n" 419 + "34.9438333 -119.6651667 671.1 \n" 420 + "34.9438333 -119.459 873.1 \n" 421 + "34.944 -119.4591667 873.2 \n" 422 + "34.9441667 -119.4101667 1064.5\n\n"); 423 424 sb 425 .append("[output file] --> the outputfile, which by default has a structure \n " 426 + "of mimicing the following pattern: \n"); 427 428 sb 429 .append(" >>> x, y, interpolated value, row level, column level <<<\n\n"); 430 431 sb 432 .append("[weight coefficient] --> the weight to be used for the \n" 433 + "interpolation, the smaller this value (or closer to 1) the \n" 434 + "smoother the grid, the larger the value the better the fit \n" 435 + "to the data. Typically, for geological data a value of 2 \n" 436 + "is used and for biological data 1.2 is commonly used.\n"); 437 438 sb.append("[search radius] --> the search radius "); 439 440 sb.append("\n\n\nNOTES\n"); 441 sb 442 .append("** if the output file extension is either .asc or .ASC an ESRI \n " 443 + "ASCII raster will be written instead of the default grid structure"); 444 445 sb.append("\n\n\nEXAMPLES\n"); 446 447 sb 448 .append("java invdist_power_isosearch2d 34.90033 35 -119.9958 -119.0037 .01 .01 -999 ../data/latLongElevation.txt.dat out.dat 1 100"); 449 sb 450 .append("java invdist_power_isosearch2d 34.90033 35 -119.9958 -119.0037 .01 .01 -999 ../data/latLongElevation.txt.dat out.asc 1 100"); 451 452 sb.append("\n\n\nFEEDBACK\n"); 453 sb.append("John Harris -- harris@numericsolutions.com"); 454 455 return sb.toString(); 456 } 457 458 public static void main(String[] args) { 459 460 try { 461 462 if (args.length != 11) { 463 System.out.println(getUsage()); 464 } else { 465 466 float xmin = (float) Double.parseDouble(args[0]); 467 float xmax = (float) Double.parseDouble(args[1]); 468 float ymin = (float) Double.parseDouble(args[2]); 469 float ymax = (float) Double.parseDouble(args[3]); 470 471 float dx = (float) Double.parseDouble(args[4]); 472 float dy = (float) Double.parseDouble(args[5]); 473 float nullval = (float) Double.parseDouble(args[6]); 474 475 String datafile = args[7]; 476 String outFile = args[8]; 477 478 float weight = (float) Double.parseDouble(args[9]); 479 float searchRadius = (float) Double.parseDouble(args[10]); 480 481 invdist_power_isosearch2d gridder = new invdist_power_isosearch2d(); 482 483 Vector template = gridder.getTemplateGrid(xmin, xmax, ymin, 484 ymax, dx, dy, nullval); 485 Vector outGrid = gridder.getInverseDistance2dGrid(template, 486 datafile, outFile, weight, searchRadius); 487 488 if (outFile.endsWith(".asc") || outFile.endsWith(".ASC")) { 489 StringBuffer out = gridder.writeArcGrid(outGrid, xmin, 490 xmax, ymin, ymax, gridder.NX, gridder.NY, dx, 491 nullval); 492 493 PrintStream outfile = new PrintStream(new FileOutputStream( 494 outFile)); 495 outfile.print(out.toString()); 496 497 } 498 } 499 500 } catch (Exception e) { 501 e.printStackTrace(); 502 } 503 504 } 505 506 // execute(dataFile, gridFile, xmin, xmax, ymin, ymax, dx, dy, nullval, 507 // weight, searchradius); 508 509 public void execute(String datafile, String outFile, float xmin, 510 float xmax, float ymin, float ymax, float dx, float dy, 511 float nullval, float weight, float searchRadius) { 512 513 try { 514 515 Vector template = getTemplateGrid(xmin, xmax, ymin, ymax, dx, dy, 516 nullval); 517 Vector outGrid = getInverseDistance2dGrid(template, datafile, 518 outFile, weight, searchRadius); 519 520 if (outFile.endsWith(".asc") || outFile.endsWith(".ASC")) { 521 StringBuffer out = writeArcGrid(outGrid, xmin, xmax, ymin, 522 ymax, NX, NY, dx, nullval); 523 524 PrintStream outfile = new PrintStream(new FileOutputStream( 525 outFile)); 526 outfile.print(out.toString()); 527 outfile.close(); // closing the file after finished. 528 } 529 530 } catch (Exception e) { 531 e.printStackTrace(); 532 } 533 534 } 535 536}