/alaspatial/src/main/java/org/ala/spatial/util/SimpleRegion.java
Java | 1665 lines | 1352 code | 58 blank | 255 comment | 171 complexity | 4ebb7388708c9748ec323c8d96270f30 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
1package org.ala.spatial.util; 2 3import java.awt.image.BufferedImage; 4import java.io.File; 5import java.io.Serializable; 6import java.util.ArrayList; 7import java.util.HashMap; 8import javax.imageio.ImageIO; 9import org.ala.spatial.analysis.cluster.SpatialCluster3; 10 11/** 12 * SimpleRegion enables point to shape intersections, where the shape 13 * is stored within SimpleRegion as a circle, bounding box or polygon. 14 * 15 * Other utilities include shape presence on a defined grid; 16 * fully present, partially present, absent. 17 * 18 * @author Adam Collins 19 */ 20public class SimpleRegion extends Object implements Serializable { 21 22 static final long serialVersionUID = -5509351896749940566L; 23 /** 24 * shape type not declared 25 */ 26 public static final int UNDEFINED = 0; 27 /** 28 * shape type bounding box; upper, lower, left, right 29 */ 30 public static final int BOUNDING_BOX = 1; 31 /** 32 * shape type circle; point and radius 33 */ 34 public static final int CIRCLE = 2; 35 /** 36 * shape type polygon; list of points as longitude, latitude pairs 37 * last point == first point 38 */ 39 public static final int POLYGON = 3; 40 /** 41 * UNDEFINED state for grid intersection output 42 * 43 * can be considered ABSENCE 44 */ 45 public static final int GI_UNDEFINED = 0; 46 /** 47 * PARTiALLy PRESENT state for grid intersection output 48 */ 49 public static final int GI_PARTIALLY_PRESENT = 1; 50 /** 51 * FULLY PRESENT state for grid intersection output 52 */ 53 public static final int GI_FULLY_PRESENT = 2; 54 /** 55 * ABSENCE state for grid intersection output 56 */ 57 public static final int GI_ABSENCE = 0; 58 /** 59 * assigned shape type 60 */ 61 int type; 62 /** 63 * points store 64 * BOUNDING_BOX = double [2][2] 65 * CIRCLE = double [1][2] 66 * POLYGON, n points (start = end) = double[n][2] 67 */ 68 float[][] points; 69 /** 70 * for point/grid to polygon intersection method 71 * 72 * polygon edges as lines of the form <code>y = a*x + b</code> 73 * lines = double [n][2] 74 * where 75 * n is number of edges 76 * value at [0] is <code>a</code> 77 * value at [1] is <code>b</code> 78 */ 79 float[][] lines2; 80 /** 81 * bounding box for types BOUNDING_BOX and POLYGON 82 * 83 * bounding_box = double [2][2] 84 * where 85 * [0][0] = minimum longitude 86 * [0][1] = minimum latitude 87 * [1][0] = maximum longitude 88 * [1][1] = maximum latitude 89 */ 90 double[][] bounding_box; //for polygons 91 /** 92 * radius for type CIRCLE in m 93 * 94 */ 95 double radius; 96 /** 97 * misc attributes 98 */ 99 HashMap<String, Object> attributes; 100 101 /** 102 * Constructor for a SimpleRegion with no shape 103 */ 104 public SimpleRegion() { 105 type = UNDEFINED; 106 } 107 108 /** 109 * gets number of points for type POLYGON 110 * 111 * note: first point = last point 112 * 113 * @return number of points as int 114 */ 115 public int getNumberOfPoints() { 116 return points.length; 117 } 118 119 /** 120 * gets the bounding box for types POLYGON and BOUNDING_BOX 121 * 122 * @return bounding box as double[2][2] 123 * with [][0] longitude and [][1] latitude 124 * minimum values at [0][], maximum values at [1][0] 125 */ 126 public double[][] getBoundingBox() { 127 return bounding_box; 128 } 129 130 /** 131 * defines the SimpleRegion as type BOUNDING_BOX 132 * 133 * @param longitude1 134 * @param latitude1 135 * @param longitude2 136 * @param latitude2 137 */ 138 public void setBox(double longitude1, double latitude1, double longitude2, double latitude2) { 139 type = BOUNDING_BOX; 140 points = new float[2][2]; 141 points[0][0] = (float) Math.min(longitude1, longitude2); 142 points[0][1] = (float) Math.min(latitude1, latitude2); 143 points[1][0] = (float) Math.max(longitude1, longitude2); 144 points[1][1] = (float) Math.max(latitude1, latitude2); 145 146 for (int i = 0; i < points.length; i++) { 147 //fix at -180 and 180 148 if (points[i][0] < -180) { 149 points[i][0] = -180; 150 } 151 if (points[i][0] > 180) { 152 points[i][0] = 180; 153 } 154 while (points[i][1] < -180) { 155 points[i][1] = -180; 156 } 157 while (points[i][1] > 180) { 158 points[i][1] = 180; 159 } 160 } 161 162 bounding_box = new double[2][2]; 163 bounding_box[0][0] = points[0][0]; 164 bounding_box[0][1] = points[0][1]; 165 bounding_box[1][0] = points[1][0]; 166 bounding_box[1][1] = points[1][1]; 167 } 168 169 /** 170 * defines the SimpleRegion as type UNDEFINED 171 */ 172 public void setNone() { 173 type = UNDEFINED; 174 } 175 176 /** 177 * defines the SimpleRegion as type CIRCLE 178 * 179 * @param longitude 180 * @param latitude 181 * @param radius_ radius of the circle in m 182 */ 183 public void setCircle(double longitude, double latitude, double radius_) { 184 type = CIRCLE; 185 points = new float[1][2]; 186 points[0][0] = (float) longitude; 187 points[0][1] = (float) latitude; 188 radius = radius_; 189 } 190 191 /** 192 * defines the SimpleRegion as type POLYGON 193 * 194 * @param points_ array of points as longitude and latiude 195 * in double [n][2] where n is the number of points 196 */ 197 public void setPolygon(double[][] points_) { 198 if (points_ != null && points_.length > 1) { 199 type = POLYGON; 200 int i; 201 202 for (i = 0; i < points_.length; i++) { 203 //fix at -180 and 180 204 if (points_[i][0] < -180) { 205 points_[i][0] = -180; 206 } 207 if (points_[i][0] > 180) { 208 points_[i][0] = 180; 209 } 210 while (points_[i][1] < -180) { 211 points_[i][1] = -180; 212 } 213 while (points_[i][1] > 180) { 214 points_[i][1] = 180; 215 } 216 } 217 218 /* copy and ensure last point == first point */ 219 int len = points_.length - 1; 220 if (points_[0][0] != points_[len][0] || points_[0][1] != points_[len][1]) { 221 points = new float[points_.length + 1][2]; 222 for (i = 0; i < points_.length; i++) { 223 points[i][0] = (float) points_[i][0]; 224 points[i][1] = (float) points_[i][1]; 225 } 226 points[points_.length][0] = (float) points_[0][0]; 227 points[points_.length][1] = (float) points_[0][1]; 228 } else { 229 points = new float[points_.length][2]; 230 for (i = 0; i < points_.length; i++) { 231 points[i][0] = (float) points_[i][0]; 232 points[i][1] = (float) points_[i][1]; 233 } 234 } 235 236 /* bounding box setup */ 237 bounding_box = new double[2][2]; 238 bounding_box[0][0] = points[0][0]; 239 bounding_box[0][1] = points[0][1]; 240 bounding_box[1][0] = points[0][0]; 241 bounding_box[1][1] = points[0][1]; 242 for (i = 1; i < points.length; i++) { 243 if (bounding_box[0][0] > points[i][0]) { 244 bounding_box[0][0] = points[i][0]; 245 } 246 if (bounding_box[1][0] < points[i][0]) { 247 bounding_box[1][0] = points[i][0]; 248 } 249 if (bounding_box[0][1] > points[i][1]) { 250 bounding_box[0][1] = points[i][1]; 251 } 252 if (bounding_box[1][1] < points[i][1]) { 253 bounding_box[1][1] = points[i][1]; 254 } 255 } 256 257 // intersection method precalculated data 258 lines2 = new float[points.length][2]; //lines[0][] is not used 259 for (i = 0; i < points.length - 1; i++) { 260 lines2[i + 1][0] = (points[i][1] - points[i + 1][1]) 261 / (points[i][0] - points[i + 1][0]); //slope 262 lines2[i + 1][1] = points[i][1] - lines2[i + 1][0] * points[i][0]; //intercept 263 } 264 } 265 } 266 267 /** 268 * defines the SimpleRegion as type POLYGON 269 * 270 * @param points_ array of points as longitude and latiude 271 * in double [n][2] where n is the number of points 272 */ 273 public void setPolygon(float[][] points_) { 274 if (points_ != null && points_.length > 1) { 275 type = POLYGON; 276 int i; 277 278 for (i = 0; i < points_.length; i++) { 279 //fix at -180 and 180 280 if (points_[i][0] < -180) { 281 points_[i][0] = -180; 282 } 283 if (points_[i][0] > 180) { 284 points_[i][0] = 180; 285 } 286 while (points_[i][1] < -180) { 287 points_[i][1] = -180; 288 } 289 while (points_[i][1] > 180) { 290 points_[i][1] = 180; 291 } 292 } 293 294 /* copy and ensure last point == first point */ 295 int len = points_.length - 1; 296 if (points_[0][0] != points_[len][0] || points_[0][1] != points_[len][1]) { 297 points = new float[points_.length + 1][2]; 298 for (i = 0; i < points_.length; i++) { 299 points[i][0] = (float) points_[i][0]; 300 points[i][1] = (float) points_[i][1]; 301 } 302 points[points_.length][0] = (float) points_[0][0]; 303 points[points_.length][1] = (float) points_[0][1]; 304 } else { 305 points = new float[points_.length][2]; 306 for (i = 0; i < points_.length; i++) { 307 points[i][0] = (float) points_[i][0]; 308 points[i][1] = (float) points_[i][1]; 309 } 310 } 311 312 /* bounding box setup */ 313 bounding_box = new double[2][2]; 314 bounding_box[0][0] = points[0][0]; 315 bounding_box[0][1] = points[0][1]; 316 bounding_box[1][0] = points[0][0]; 317 bounding_box[1][1] = points[0][1]; 318 for (i = 1; i < points.length; i++) { 319 if (bounding_box[0][0] > points[i][0]) { 320 bounding_box[0][0] = points[i][0]; 321 } 322 if (bounding_box[1][0] < points[i][0]) { 323 bounding_box[1][0] = points[i][0]; 324 } 325 if (bounding_box[0][1] > points[i][1]) { 326 bounding_box[0][1] = points[i][1]; 327 } 328 if (bounding_box[1][1] < points[i][1]) { 329 bounding_box[1][1] = points[i][1]; 330 } 331 } 332 333 // intersection method precalculated data 334 lines2 = new float[points.length][2]; //lines[0][] is not used 335 for (i = 0; i < points.length - 1; i++) { 336 lines2[i + 1][0] = (points[i][1] - points[i + 1][1]) 337 / (points[i][0] - points[i + 1][0]); //slope 338 lines2[i + 1][1] = points[i][1] - lines2[i + 1][0] * points[i][0]; //intercept 339 } 340 } 341 } 342 343 /** 344 * gets points of a polygon only 345 * 346 * @return points of this object if it is a polygon as double[][] 347 * otherwise returns null. 348 */ 349 public float[][] getPoints() { 350 return points; 351 } 352 353 /** 354 * returns true when the point provided is within the SimpleRegion 355 * 356 * note: type UNDEFINED implies no boundary, always returns true. 357 * 358 * @param longitude 359 * @param latitude 360 * @return true iff point is within or on the edge of this SimpleRegion 361 */ 362 public boolean isWithin(double longitude, double latitude) { 363 switch (type) { 364 case 0: 365 /* no region defined, must be within this absence of a boundary */ 366 return true; 367 case 1: 368 /* return for bounding box */ 369 return (longitude <= points[1][0] && longitude >= points[0][0] 370 && latitude <= points[1][1] && latitude >= points[0][1]); 371 case 2: 372 /* TODO: fix to use radius units m not degrees */ 373 double x = longitude - points[0][0]; 374 double y = latitude - points[0][1]; 375 return Math.sqrt(x * x + y * y) <= radius; 376 case 3: 377 /* determine for Polygon */ 378 return isWithinPolygon(longitude, latitude); 379 } 380 return false; 381 } 382 383 public boolean isWithin_EPSG900913(double longitude, double latitude) { 384 switch (type) { 385 case 0: 386 /* no region defined, must be within this absence of a boundary */ 387 return true; 388 case 1: 389 /* return for bounding box */ 390 return (longitude <= points[1][0] && longitude >= points[0][0] 391 && latitude <= points[1][1] && latitude >= points[0][1]); 392 case 2: 393 /* TODO: fix to use radius units m not degrees */ 394 double x = longitude - points[0][0]; 395 double y = latitude - points[0][1]; 396 return Math.sqrt(x * x + y * y) <= radius; 397 case 3: 398 /* determine for Polygon */ 399 return isWithinPolygon_EPSG900913(longitude, latitude); 400 } 401 return false; 402 } 403 404 /** 405 * returns true when point is within the polygon 406 * 407 * method: 408 * treat as segments with target long in the middle: 409 * 410 * 411 * __-1__|___1_ 412 * | 413 * 414 * 415 * iterate through points and count number of latitude axis crossing where 416 * crossing is > latitude. 417 * 418 * point is inside of area when number of crossings is odd; 419 * 420 * point is on a polygon edge return true 421 * 422 * @param longitude 423 * @param latitude 424 * @return true iff longitude and latitude point is on edge or within polygon 425 */ 426 private boolean isWithinPolygon(double longitude, double latitude) { 427 // bounding box test 428 if (longitude <= bounding_box[1][0] && longitude >= bounding_box[0][0] 429 && latitude <= bounding_box[1][1] && latitude >= bounding_box[0][1]) { 430 431 //initial segment 432 boolean segment = points[0][0] > longitude; 433 434 double y; 435 int i; 436 int len = points.length; 437 int score = 0; 438 439 for (i = 1; i < len; i++) { 440 // is it in a new segment? 441 if ((points[i][0] > longitude) != segment) { 442 //lat value at line crossing > target point 443 y = lines2[i][0] * longitude + lines2[i][1]; 444 if (y > latitude) { 445 score++; 446 } else if (y == latitude) { 447 //line crossing 448 return true; 449 } 450 451 segment = !segment; 452 } else if (points[i][0] == longitude && points[i][1] == latitude) { 453 //point on point 454 return true; 455 } 456 } 457 return (score % 2 != 0); 458 } 459 return false; //not within bounding box 460 } 461 462 private boolean isWithinPolygon_EPSG900913(double longitude, double latitude) { 463 SpatialCluster3 sc = new SpatialCluster3(); 464 465 // bounding box test 466 if (longitude <= bounding_box[1][0] && longitude >= bounding_box[0][0] 467 && latitude <= bounding_box[1][1] && latitude >= bounding_box[0][1]) { 468 469 //initial segment 470 int longitudePx = sc.convertLngToPixel(longitude); 471 boolean segment = sc.convertLngToPixel(points[0][0]) > longitudePx; 472 473 int y; 474 int i; 475 int len = points.length; 476 int score = 0; 477 478 for (i = 1; i < len; i++) { 479 // is it in a new segment? 480 if ((sc.convertLngToPixel(points[i][0]) > longitudePx) != segment) { 481 //lat value at line crossing > target point 482 y = (int)((longitudePx - sc.convertLngToPixel(points[i][0])) 483 * ((sc.convertLatToPixel(points[i][1]) - sc.convertLatToPixel(points[i - 1][1])) 484 / (double) (sc.convertLngToPixel(points[i][0]) - sc.convertLngToPixel(points[i - 1][0]))) 485 + sc.convertLatToPixel(points[i][1])); 486 if (y > sc.convertLatToPixel(latitude)) { 487 score++; 488 } else if (y == sc.convertLatToPixel(latitude)) { 489 //line crossing 490 return true; 491 } 492 493 segment = !segment; 494 } else if (points[i][0] == longitude && points[i][1] == latitude) { 495 //point on point 496 return true; 497 } 498 } 499 return (score % 2 != 0); 500 } 501 return false; //not within bounding box 502 } 503 504 /** 505 * determines overlap with a grid 506 * 507 * for type POLYGON 508 * when <code>three_state_map</code> is not null populate it with one of: 509 * GI_UNDEFINED 510 * GI_PARTIALLY_PRESENT 511 * GI_FULLY_PRESENT 512 * GI_ABSENCE 513 * 514 * @param longitude1 515 * @param latitude1 516 * @param longitude2 517 * @param latitude2 518 * @param xres number of longitude segements as int 519 * @param yres number of latitude segments as int 520 * @return (x,y) as double [][2] for each grid cell at least partially falling 521 * within the specified region of the specified resolution beginning at 0,0 522 * for minimum longitude and latitude through to xres,yres for maximums 523 */ 524 public int[][] getOverlapGridCells(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) { 525 int[][] cells = null; 526 switch (type) { 527 case 0: 528 break; 529 case 1: 530 cells = getOverlapGridCells_Box(longitude1, latitude1, longitude2, latitude2, width, height, bounding_box, three_state_map, noCellsReturned); 531 break; 532 case 2: 533 break; /* TODO: circle grid */ 534 case 3: 535 cells = getOverlapGridCells_Polygon(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map, noCellsReturned); 536 } 537 538 return cells; 539 } 540 541 /** 542 * stacks PARTIALLY_PRESENT shape outline onto three_state_map 543 * 544 * @param longitude1 545 * @param latitude1 546 * @param longitude2 547 * @param latitude2 548 * @param width 549 * @param height 550 * @param three_state_map 551 * @param noCellsReturned 552 */ 553 public void getOverlapGridCells_Acc(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map) { 554 switch (type) { 555 case 0: 556 break; 557 case 1: 558 getOverlapGridCells_Box_Acc(longitude1, latitude1, longitude2, latitude2, width, height, bounding_box, three_state_map); 559 break; 560 case 2: 561 break; /* TODO: circle grid */ 562 case 3: 563 getOverlapGridCells_Polygon_Acc(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map); 564 } 565 } 566 567 public int[][] getOverlapGridCells(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map) { 568 return getOverlapGridCells(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map, false); 569 } 570 571 /** 572 * determines overlap with a grid for a BOUNDING_BOX 573 * 574 * @param longitude1 575 * @param latitude1 576 * @param longitude2 577 * @param latitude2 578 * @param xres number of longitude segements as int 579 * @param yres number of latitude segments as int 580 * @param bb bounding box as double[2][2] with [][0] as longitude, [][1] as latitude, 581 * [0][] as minimum values, [1][] as maximum values 582 * @return (x,y) as double [][2] for each grid cell at least partially falling 583 * within the specified region of the specified resolution beginning at 0,0 584 * for minimum longitude and latitude through to xres,yres for maximums 585 */ 586 public int[][] getOverlapGridCells_Box(double longitude1, double latitude1, 587 double longitude2, double latitude2, int width, int height, double[][] bb, byte[][] three_state_map, boolean noCellsReturned) { 588 589 double xstep = Math.abs(longitude2 - longitude1) / (double) width; 590 double ystep = Math.abs(latitude2 - latitude1) / (double) height; 591 592 //double maxlong = Math.max(longitude1, longitude2); 593 double minlong = Math.min(longitude1, longitude2); 594 //double maxlat = Math.max(latitude1, latitude2); 595 double minlat = Math.min(latitude1, latitude2); 596 597 //setup minimums from bounding box (TODO: should this have -1 on steps?) 598 int xstart = (int) Math.floor((bb[0][0] - minlong) / xstep); 599 int ystart = (int) Math.floor((bb[0][1] - minlat) / ystep); 600 int xend = (int) Math.ceil((bb[1][0] - minlong) / xstep); 601 int yend = (int) Math.ceil((bb[1][1] - minlat) / ystep); 602 if (xstart < 0) { 603 xstart = 0; 604 } 605 if (ystart < 0) { 606 ystart = 0; 607 } 608 if (xend > width) { 609 xend = width; 610 } 611 if (yend > height) { 612 yend = height; 613 } 614 615 // fill data with cell coordinates 616 int out_width = xend - xstart; 617 int out_height = yend - ystart; 618 int j, i, p = 0; 619 int[][] data = null; 620 if (!noCellsReturned) { 621 data = new int[out_width * out_height][2]; 622 if (three_state_map == null) { 623 for (j = ystart; j < yend; j++) { 624 for (i = xstart; i < xend; i++) { 625 data[p][0] = i; 626 data[p][1] = j; 627 p++; 628 } 629 } 630 } else { 631 for (j = ystart; j < yend; j++) { 632 for (i = xstart; i < xend; i++) { 633 data[p][0] = i; 634 data[p][1] = j; 635 three_state_map[j][i] = SimpleRegion.GI_FULLY_PRESENT; 636 p++; 637 } 638 } 639 //set three state map edges to partially present 640 if (xstart < xend && xend > 0) { 641 for (j = ystart; j < yend; j++) { 642 three_state_map[j][xstart] = SimpleRegion.GI_PARTIALLY_PRESENT; 643 three_state_map[j][xend - 1] = SimpleRegion.GI_PARTIALLY_PRESENT; 644 } 645 } 646 if (ystart < yend && yend > 0) { 647 for (i = xstart; i < xend; i++) { 648 three_state_map[ystart][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 649 three_state_map[yend - 1][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 650 } 651 } 652 653 //no need to set SimpleRegion.GI_ABSENCE 654 } 655 } else { 656 if (three_state_map == null) { 657 for (j = ystart; j < yend; j++) { 658 for (i = xstart; i < xend; i++) { 659 data[p][0] = i; 660 data[p][1] = j; 661 p++; 662 } 663 } 664 } else { 665 for (j = ystart; j < yend; j++) { 666 for (i = xstart; i < xend; i++) { 667 three_state_map[j][i] = SimpleRegion.GI_FULLY_PRESENT; 668 } 669 } 670 //set three state map edges to partially present 671 if (xstart < xend && xend > 0) { 672 for (j = ystart; j < yend; j++) { 673 three_state_map[j][xstart] = SimpleRegion.GI_PARTIALLY_PRESENT; 674 three_state_map[j][xend - 1] = SimpleRegion.GI_PARTIALLY_PRESENT; 675 } 676 } 677 if (ystart < yend && yend > 0) { 678 for (i = xstart; i < xend; i++) { 679 three_state_map[ystart][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 680 three_state_map[yend - 1][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 681 } 682 } 683 } 684 } 685 return data; 686 } 687 688 /** 689 * defines a region by a points string, POLYGON only 690 * 691 * TODO: define better format for parsing, including BOUNDING_BOX and CIRCLE 692 * 693 * @param pointsString points separated by ',' with longitude and latitude separated by ':' 694 * @return SimpleRegion object 695 */ 696 public static SimpleRegion parseSimpleRegion(String pointsString) { 697 if (pointsString.equalsIgnoreCase("none")) { 698 return null; 699 } 700 SimpleRegion simpleregion = new SimpleRegion(); 701 ArrayList<Double> points = new ArrayList<Double>(); 702 703 int pos; 704 int lastpos = 0; 705 while((pos = Math.min(pointsString.indexOf(',', lastpos), pointsString.indexOf(' ', lastpos))) > 0) { 706 try { 707 points.add(Double.parseDouble(pointsString.substring(lastpos, pos))); 708 } catch (Exception e) { 709 points.add(0.0); 710 } 711 lastpos = pos + 1; 712 } 713 //one coordinate pair left 714 pos = pointsString.indexOf(' ', lastpos); 715 try { 716 points.add(Double.parseDouble(pointsString.substring(lastpos, pos))); 717 lastpos = pos + 1; 718 } catch (Exception e) { 719 points.add(0.0); 720 } 721 try { 722 points.add(Double.parseDouble(pointsString.substring(lastpos, pointsString.length()))); 723 } catch (Exception e) { 724 points.add(0.0); 725 } 726 727 //test for box 728 // get min/max long/lat 729 // each point has only one identical lat or long to previous point 730 // 4 or 5 points (start and end points may be identical) 731 if (((points.size() == 8 && (points.get(0) != points.get(6) || points.get(1) != points.get(7))) 732 || (points.size() == 5 && points.get(0) == points.get(8) 733 && points.get(1) == points.get(9)))) { 734 735 //get min/max long/lat 736 double minlong = 0, minlat = 0, maxlong = 0, maxlat = 0; 737 for (int i = 0; i < points.size(); i+=2) { 738 if (i == 0 || minlong > points.get(i)) { 739 minlong = points.get(i); 740 } 741 if (i == 0 || maxlong < points.get(i)) { 742 maxlong = points.get(i); 743 } 744 if (i == 0 || minlat > points.get(i+1)) { 745 minlat = points.get(i+1); 746 } 747 if (i == 0 || maxlat < points.get(i+1)) { 748 maxlat = points.get(i+1); 749 } 750 } 751 752 // each point has only one identical lat or long to previous point 753 int prev_idx = 6; 754 int i = 0; 755 for (i = 0; i < 8; i+=2) { 756 if ((points.get(i) == points.get(prev_idx)) 757 == (points.get(i+1) == points.get(prev_idx+1))) { 758 break; 759 } 760 prev_idx = i; 761 } 762 //it is a box if no 'break' occurred 763 if (i == 8) { 764 simpleregion.setBox(minlong, minlat, maxlong, maxlat); 765 return simpleregion; 766 } 767 } 768 769 double [][] pointsArray = new double[points.size()/2][2]; 770 for(int i=0;i<points.size();i+=2) { 771 pointsArray[i/2][0] = points.get(i); 772 pointsArray[i/2][1] = points.get(i+1); 773 } 774 simpleregion.setPolygon(pointsArray); 775 return simpleregion; 776 } 777 778 public double getWidth() { 779 return bounding_box[1][0] - bounding_box[0][0]; 780 } 781 782 public double getHeight() { 783 return bounding_box[1][1] - bounding_box[0][1]; 784 } 785 786 public int getType() { 787 return type; 788 } 789 790 public void setAttribute(String name, Object value) { 791 if (attributes == null) { 792 attributes = new HashMap<String, Object>(); 793 } 794 attributes.put(name, value); 795 } 796 797 public Object getAttribute(String name) { 798 if (attributes != null) { 799 return attributes.get(name); 800 } 801 return null; 802 } 803 804 public void saveGridAsImage(byte[][] three_state_map) { 805 try { 806 long t1 = System.currentTimeMillis(); 807 BufferedImage bi = new BufferedImage(three_state_map[0].length, three_state_map.length, BufferedImage.TYPE_INT_RGB); 808 for (int i = 0; i < three_state_map.length; i++) { 809 for (int j = 0; j < three_state_map[i].length; j++) { 810 if (three_state_map[i][j] == 0) { 811 bi.setRGB(j, (three_state_map.length - 1 - i), 0xFFFFFF); 812 } else if (three_state_map[i][j] == 1) { 813 bi.setRGB(j, (three_state_map.length - 1 - i), 0x99FF99); 814 } else if (three_state_map[i][j] == 2) { 815 bi.setRGB(j, (three_state_map.length - 1 - i), 0x9999FF); 816 } 817 } 818 } 819 820 ImageIO.write(bi, "png", File.createTempFile("grd", ".png", new File("d:\\"))); 821 System.out.println("save grid in " + (System.currentTimeMillis() - t1) + "ms"); 822 823 } catch (Exception e) { 824 e.printStackTrace(); 825 } 826 } 827 828 /** 829 * determines overlap with a grid for POLYGON 830 * 831 * when <code>three_state_map</code> is not null populate it with one of: 832 * GI_UNDEFINED 833 * GI_PARTIALLY_PRESENT 834 * GI_FULLY_PRESENT 835 * GI_ABSENCE 836 * 837 * 1. Get 3state mask and fill edge passes as 'partial'. 838 * then 839 * 3. Test 0,0 then progress across vert raster until finding cells[][] entry 840 * 4. Repeat from (3). 841 * 842 * @param longitude1 843 * @param latitude1 844 * @param longitude2 845 * @param latitude2 846 * @param xres number of longitude segements as int 847 * @param yres number of latitude segments as int 848 * @return (x,y) as double [][2] for each grid cell at least partially falling 849 * within the specified region of the specified resolution beginning at 0,0 850 * for minimum longitude and latitude through to xres,yres for maximums 851 */ 852 public int[][] getOverlapGridCells_Polygon(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) { 853 int i, j; 854 if (three_state_map == null) { 855 three_state_map = new byte[height][width]; 856 } 857 858 double divx = (longitude2 - longitude1) / width; 859 double divy = (latitude2 - latitude1) / height; 860 861 //to cells 862 int x, y, xend, yend, xDirection, icross; 863 double xcross, endlat, dx1, dx2, dy1, dy2, slope, intercept; 864 for (j = 1; j < points.length; j++) { 865 if (points[j][1] < points[j - 1][1]) { 866 dx1 = points[j][0]; 867 dy1 = points[j][1]; 868 dx2 = points[j - 1][0]; 869 dy2 = points[j - 1][1]; 870 } else { 871 dx2 = points[j][0]; 872 dy2 = points[j][1]; 873 dx1 = points[j - 1][0]; 874 dy1 = points[j - 1][1]; 875 } 876 x = (int) ((dx1 - longitude1) / divx); 877 y = (int) ((dy1 - latitude1) / divy); 878 xend = (int) ((dx2 - longitude1) / divx); 879 yend = (int) ((dy2 - latitude1) / divy); 880 881 if (y >= 0 && y < height && x >= 0 && x < width) { 882 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 883 } 884 885 if (x == xend && y == yend) { 886 continue; 887 } 888 889 xDirection = (x < xend) ? 1 : -1; 890 891 slope = (dy1 - dy2) / (dx1 - dx2); 892 intercept = dy1 - slope * dx1; 893 894 if (x == xend) { 895 //vertical line 896 while (y != yend) { 897 y++; 898 if (y >= 0 && y < height && x >= 0 && x < width) { 899 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 900 } 901 } 902 } else if (y == yend) { 903 //horizontal line 904 while (x != xend) { 905 x += xDirection; 906 if (y >= 0 && y < height && x >= 0 && x < width) { 907 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 908 } 909 } 910 } else { //sloped line 911 endlat = dy2; 912 for (double k = (y + 1) * divy + latitude1; k < endlat; k += divy) { 913 //move in yDirection to get x 914 xcross = (k - intercept) / slope; 915 icross = (int) ((xcross - longitude1) / divx); 916 917 while (x != icross && x != xend) { 918 x += xDirection; 919 if (y >= 0 && y < height && x >= 0 && x < width) { 920 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 921 } 922 } 923 924 if (y != yend) { 925 y++; 926 if (y >= 0 && y < height && x >= 0 && x < width) { 927 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 928 } 929 } 930 } 931 932 //finish horizontal line 933 while (x != xend) { 934 x += xDirection; 935 if (y >= 0 && y < height && x >= 0 && x < width) { 936 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 937 } 938 } 939 } 940 } 941 942 //do raster check 943 int[][] data = new int[width * height][2]; 944 boolean cellsReturned = !noCellsReturned; 945 int p = 0; 946 for (j = 0; j < three_state_map[0].length; j++) { 947 for (i = 0; i < three_state_map.length; i++) { 948 if (three_state_map[i][j] == GI_PARTIALLY_PRESENT) { 949 //if it is partially present, do nothing 950 } else if ((j == 0 || three_state_map[i][j - 1] == GI_PARTIALLY_PRESENT)) { 951 if (i > 0 952 && (three_state_map[i - 1][j] == GI_FULLY_PRESENT 953 || three_state_map[i - 1][j] == GI_ABSENCE)) { 954 //use same as LHS 955 three_state_map[i][j] = three_state_map[i - 1][j]; 956 } else if (isWithin(j * divx + divx / 2 + longitude1, i * divy + divy / 2 + latitude1)) { 957 //if the previous was partially present, test 958 three_state_map[i][j] = GI_FULLY_PRESENT; 959 } //else absent 960 } else { 961 //if the previous was fully present, repeat 962 //if the previous was absent, repeat 963 three_state_map[i][j] = three_state_map[i][j - 1]; 964 } 965 966 //apply to cells; 967 if (cellsReturned && three_state_map[i][j] != GI_UNDEFINED) { //undefined == absence 968 data[p][0] = j; 969 data[p][1] = i; 970 p++; 971 } 972 } 973 } 974 return java.util.Arrays.copyOfRange(data, 0, p); 975 } 976 977 public void getOverlapGridCells_Polygon_Acc(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map) { 978 int i, j; 979 980 double divx = (longitude2 - longitude1) / width; 981 double divy = (latitude2 - latitude1) / height; 982 983 //to cells 984 int x, y, xend, yend, xDirection, icross; 985 double xcross, endlat, dx1, dx2, dy1, dy2, slope, intercept; 986 for (j = 1; j < points.length; j++) { 987 if (points[j][1] < points[j - 1][1]) { 988 dx1 = points[j][0]; 989 dy1 = points[j][1]; 990 dx2 = points[j - 1][0]; 991 dy2 = points[j - 1][1]; 992 } else { 993 dx2 = points[j][0]; 994 dy2 = points[j][1]; 995 dx1 = points[j - 1][0]; 996 dy1 = points[j - 1][1]; 997 } 998 x = (int) ((dx1 - longitude1) / divx); 999 y = (int) ((dy1 - latitude1) / divy); 1000 xend = (int) ((dx2 - longitude1) / divx); 1001 yend = (int) ((dy2 - latitude1) / divy); 1002 1003 if (x == xend && y == yend) { 1004 continue; 1005 } 1006 1007 xDirection = (x < xend) ? 1 : -1; 1008 1009 slope = (dy1 - dy2) / (dx1 - dx2); 1010 intercept = dy1 - slope * dx1; 1011 1012 if (y >= 0 && y < height && x >= 0 && x < width) { 1013 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1014 } 1015 1016 if (x == xend) { 1017 //vertical line 1018 while (y != yend) { 1019 y++; 1020 if (y >= 0 && y < height && x >= 0 && x < width) { 1021 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1022 } 1023 } 1024 } else if (y == yend) { 1025 //horizontal line 1026 while (x != xend) { 1027 x += xDirection; 1028 if (y >= 0 && y < height && x >= 0 && x < width) { 1029 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1030 } 1031 } 1032 } else { //sloped line 1033 endlat = dy2; 1034 for (double k = (y + 1) * divy + latitude1; k < endlat; k += divy) { 1035 //move in yDirection to get x 1036 xcross = (k - intercept) / slope; 1037 icross = (int) ((xcross - longitude1) / divx); 1038 1039 while (x != icross && x != xend) { 1040 x += xDirection; 1041 if (y >= 0 && y < height && x >= 0 && x < width) { 1042 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1043 } 1044 } 1045 1046 if (y != yend) { 1047 y++; 1048 if (y >= 0 && y < height && x >= 0 && x < width) { 1049 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1050 } 1051 } 1052 } 1053 1054 //finish horizontal line 1055 while (x != xend) { 1056 x += xDirection; 1057 if (y >= 0 && y < height && x >= 0 && x < width) { 1058 three_state_map[y][x] = GI_PARTIALLY_PRESENT; 1059 } 1060 } 1061 } 1062 } 1063 } 1064 1065 public int[][] fillAccMask(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) { 1066 double divx = (longitude2 - longitude1) / width; 1067 double divy = (latitude2 - latitude1) / height; 1068 1069 int i, j; 1070 //do raster check 1071 int[][] data = null; 1072 boolean cellsReturned = !noCellsReturned; 1073 if (cellsReturned) { 1074 data = new int[width * height][2]; 1075 } 1076 int p = 0; 1077 for (j = 0; j < three_state_map[0].length; j++) { 1078 for (i = 0; i < three_state_map.length; i++) { 1079 if (three_state_map[i][j] == GI_PARTIALLY_PRESENT) { 1080 //if it is partially present, do nothing 1081 } else if ((j == 0 || three_state_map[i][j - 1] == GI_PARTIALLY_PRESENT)) { 1082 if (i > 0 1083 && (three_state_map[i - 1][j] == GI_FULLY_PRESENT 1084 || three_state_map[i - 1][j] == GI_ABSENCE)) { 1085 //use same as LHS 1086 three_state_map[i][j] = three_state_map[i - 1][j]; 1087 } else if (isWithin(j * divx + divx / 2 + longitude1, i * divy + divy / 2 + latitude1)) { 1088 //if the previous was partially present, test 1089 three_state_map[i][j] = GI_FULLY_PRESENT; 1090 } //else absent 1091 } else { 1092 //if the previous was fully present, repeat 1093 //if the previous was absent, repeat 1094 three_state_map[i][j] = three_state_map[i][j - 1]; 1095 } 1096 1097 //apply to cells; 1098 if (cellsReturned && three_state_map[i][j] != GI_UNDEFINED) { //undefined == absence 1099 data[p][0] = j; 1100 data[p][1] = i; 1101 p++; 1102 } 1103 } 1104 } 1105 if (data != null) { 1106 data = java.util.Arrays.copyOf(data, p); 1107 } 1108 return data; 1109 } 1110 1111 private void getOverlapGridCells_Box_Acc(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, double[][] bb, byte[][] three_state_map) { 1112 double xstep = Math.abs(longitude2 - longitude1) / (double) width; 1113 double ystep = Math.abs(latitude2 - latitude1) / (double) height; 1114 1115 //double maxlong = Math.max(longitude1, longitude2); 1116 double minlong = Math.min(longitude1, longitude2); 1117 //double maxlat = Math.max(latitude1, latitude2); 1118 double minlat = Math.min(latitude1, latitude2); 1119 1120 //setup minimums from bounding box (TODO: should this have -1 on steps?) 1121 int xstart = (int) Math.floor((bb[0][0] - minlong) / xstep); 1122 int ystart = (int) Math.floor((bb[0][1] - minlat) / ystep); 1123 int xend = (int) Math.ceil((bb[1][0] - minlong) / xstep); 1124 int yend = (int) Math.ceil((bb[1][1] - minlat) / ystep); 1125 if (xstart < 0) { 1126 xstart = 0; 1127 } 1128 if (ystart < 0) { 1129 ystart = 0; 1130 } 1131 if (xend > width) { 1132 xend = width; 1133 } 1134 if (yend > height) { 1135 yend = height; 1136 } 1137 1138 // fill data with cell coordinates 1139 //int out_width = xend - xstart; 1140 //int out_height = yend - ystart; 1141 int j, i, p = 0; 1142 //int[][] data = null; 1143 1144 //set three state map edges to partially present 1145 if (xstart < xend && xend > 0) { 1146 for (j = ystart; j < yend; j++) { 1147 three_state_map[j][xstart] = SimpleRegion.GI_PARTIALLY_PRESENT; 1148 three_state_map[j][xend - 1] = SimpleRegion.GI_PARTIALLY_PRESENT; 1149 } 1150 } 1151 if (ystart < yend && yend > 0) { 1152 for (i = xstart; i < xend; i++) { 1153 three_state_map[ystart][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 1154 three_state_map[yend - 1][i] = SimpleRegion.GI_PARTIALLY_PRESENT; 1155 } 1156 } 1157 } 1158 1159 public int[][] getOverlapGridCells_EPSG900913(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) { 1160 int[][] cells = null; 1161 switch (type) { 1162 case 0: 1163 break; 1164 case 1: 1165 cells = getOverlapGridCells_Box(longitude1, latitude1, longitude2, latitude2, width, height, bounding_box, three_state_map, noCellsReturned); 1166 break; 1167 case 2: 1168 break; /* TODO: circle grid */ 1169 case 3: 1170 cells = getOverlapGridCells_Polygon_EPSG900913(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map, noCellsReturned); 1171 } 1172 1173 return cells; 1174 } 1175 1176 /** 1177 * stacks PARTIALLY_PRESENT shape outline onto three_state_map 1178 * 1179 * @param longitude1 1180 * @param latitude1 1181 * @param longitude2 1182 * @param latitude2 1183 * @param width 1184 * @param height 1185 * @param three_state_map 1186 * @param noCellsReturned 1187 */ 1188 public void getOverlapGridCells_Acc_EPSG900913(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map) { 1189 switch (type) { 1190 case 0: 1191 break; 1192 case 1: 1193 getOverlapGridCells_Box_Acc(longitude1, latitude1, longitude2, latitude2, width, height, bounding_box, three_state_map); 1194 break; 1195 case 2: 1196 break; /* TODO: circle grid */ 1197 case 3: 1198 getOverlapGridCells_Polygon_Acc_EPSG900913(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map); 1199 } 1200 } 1201 1202 public int[][] getOverlapGridCells_EPSG900913(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map) { 1203 return getOverlapGridCells_EPSG900913(longitude1, latitude1, longitude2, latitude2, width, height, three_state_map, false); 1204 } 1205 1206 public int[][] getOverlapGridCells_Polygon_EPSG900913(double olongitude1, double olatitude1, double olongitude2, double olatitude2, int owidth, int oheight, byte[][] three_state_map, boolean noCellsReturned) { 1207 int i, j; 1208 if (three_state_map == null) { 1209 three_state_map = new byte[oheight][owidth]; 1210 } 1211 1212 SpatialCluster3 sc = new SpatialCluster3(); 1213 1214 int longitude1 = sc.convertLngToPixel(olongitude1); 1215 int longitude2 = sc.convertLngToPixel(olongitude2); 1216 int latitude1 = sc.convertLatToPixel(olatitude2); 1217 int latitude2 = sc.convertLatToPixel(olatitude1); 1218 int scale = 100; //if it is too small the 'fill' operation is messed up 1219 int width = owidth * scale; 1220 int height = oheight * scale; 1221 1222 double divx = (longitude2 - longitude1) / width; 1223 double divy = (latitude2 - latitude1) / height; 1224 double odivx = (olongitude2 - olongitude1) / owidth; 1225 double odivy = (olatitude2 - olatitude1) / oheight; 1226 1227 int oy, ox; 1228 1229 //to cells 1230 int x, y, xend, yend, xDirection, icross; 1231 double slope, intercept; 1232 int xcross, endlat, dx1, dx2, dy1, dy2; 1233 for (j = 1; j < points.length; j++) { 1234 if (points[j][1] > points[j - 1][1]) { 1235 dx1 = sc.convertLngToPixel(points[j][0]); 1236 dy1 = sc.convertLatToPixel(points[j][1]); 1237 dx2 = sc.convertLngToPixel(points[j - 1][0]); 1238 dy2 = sc.convertLatToPixel(points[j - 1][1]); 1239 } else { 1240 dx2 = sc.convertLngToPixel(points[j][0]); 1241 dy2 = sc.convertLatToPixel(points[j][1]); 1242 dx1 = sc.convertLngToPixel(points[j - 1][0]); 1243 dy1 = sc.convertLatToPixel(points[j - 1][1]); 1244 } 1245 x = (int) ((dx1 - longitude1) / divx); 1246 y = (int) ((dy1 - latitude1) / divy); 1247 xend = (int) ((dx2 - longitude1) / divx); 1248 yend = (int) ((dy2 - latitude1) / divy); 1249 1250 if (y >= 0 && y < height && x >= 0 && x < width) { 1251 oy = (int) ((sc.convertPixelToLat((int) (y * divy + latitude1)) - olatitude1) / odivy); 1252 ox = (int) ((sc.convertPixelToLng((int) (x * divx + longitude1)) - olongitude1) / odivx); 1253 if (oy >= oheight) { 1254 oy = oheight - 1; 1255 } 1256 if (ox >= owidth) { 1257 ox = owidth - 1; 1258 } 1259 if (oy < 0) { 1260 oy = 0; 1261 } 1262 if (ox < 0) { 1263 ox = 0; 1264 } 1265 three_state_map[oy][ox] = GI_PARTIALLY_PRESENT; 1266 } 1267 1268 if (x == xend && y == yend) { 1269 continue; 1270 } 1271 1272 xDirection = (x < xend) ? 1 : -1; 1273 1274 slope = (dy1 - dy2) / (double) (dx1 - dx2); 1275 intercept = (double) (dy1 - slope * dx1); 1276 1277 if (x == xend) { 1278 //vertical line 1279 while (y != yend) { 1280 y++; 1281 if (y >= 0 && y < height && x >= 0 && x < width) { 1282 oy = (int) ((sc.convertPixelToLat((int) (y * divy + latitude1)) - olatitude1) / odivy); 1283 ox = (int) ((sc.convertPixelToLng((int) (x * divx + longitude1)) - olongitude1) / odivx); 1284 if (oy >= oheight) { 1285 oy = oheight - 1; 1286 } 1287 if (ox >= owidth) { 1288 ox = owidth - 1; 1289 } 1290 if (oy < 0) { 1291 oy = 0; 1292 } 1293 if (ox < 0) { 1294 ox = 0; 1295 } 1296 three_state_map[oy][ox] = GI_PARTIALLY_PRESENT; 1297 } 1298 } 1299 } else if (y == yend) { 1300 //horizontal line 1301 while (x != xend) { 1302 x += xDirection; 1303 if (y >= 0 && y < height && x >= 0 && x < width) { 1304 oy = (int) ((sc.convertPixelToLat((int) (y * divy + latitude1)) - olatitude1) / odivy); 1305 ox = (int) ((sc.convertPixelToLng((int) (x * divx + longitude1)) - olongitude1) / odivx); 1306 if (oy >= oheight) { 1307 oy = oheight - 1; 1308 } 1309 if (ox >= owidth) { 1310 ox = owidth - 1; 1311 } 1312 if (oy < 0) { 1313 oy = 0; 1314 } 1315 if (ox < 0) { 1316 ox = 0; 1317 } 1318 …
Large files files are truncated, but you can click here to view the full file