Información espacial en formato ráster en R
Existen varios paquetes en R que permiten manejar datos espaciales, ya sea en formato de vector o raster. Algunos de ellos incluyen: sp, rgdal, rgeos, sf, stars, raster, terra. Sin embargo, en este curso nos enfocaremos en utilizar sf para el manejo de información vectorial y stars para el manejo de rasters.
Introducción a stars
stars es un paquete creado para trabajar de manera sencilla con información raster y fue creado por el mismo grupo de trabajo que creo sf, así que muchas funciones y sintaxis son compartidas con dicho paquete.
Rasters
El objeto básico en stars son los spatiotemporal arrays (stars), es decir arreglos espacio temporales. Recordemos que un arreglo corresponde a cualquier objeto organizado en dimensiones. Por ejemplo, una arreglo de dos dimensiones podría corresponder a una matriz que cuenta con una dimensión vertical (filas) y horizontal (columnas). De igual manera, un arreglo de tres dimensiones podría corresponder a una imagen raster con dos dimensiones espaciales (vertical y horizontal) y una dimensión espectral (bandas). Adicionalmente, un arreglo con cuatro dimensiones podría incluir una dimensión temporal.
Al igual que en sf cada raster va a estar asociado a un sistema de coordenadas de referencia (CRS). Por lo cual, dicho CRS indicará la proyección de los datos y el datum.
Estructura de datos stars
Los datos stars van a estar organizados en un formato de array, es decir, de una tabla o cuadro de datos de n dimensiones. Para ver un ejemplos carguemos un archivo raster.
library(stars)## Loading required package: abind## Loading required package: sf## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUEtif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
x## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## L7_ETMs.tif 1 54 69 68.91242 86 255
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [x]
## y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [y]
## band 1 6 NA NA NA NA NULLAnalicemos qué dice esta información. Primero nos dice que se trata de un objeto stars de tres dimensiones y un atributo. Además, nos da algunas estadísticas de la información que contiene. A continuación, podemos ver las tres dimensiones: x, y y band, es decir, dos dimensiones espaciales y una espectral. Adicionalmente, se puede consultar el número de celdas en cada dimensión (from y to), la coordenada inicial en cada dimensión (offset), el tamaño de la celda en cada dimensión (delta), el sistema de referencia en el que se encuentra la información (refsys), si corresponde a puntos (point) y los valores de secuencia en caso de que la dimensión no sea regular (e.g., geometrías).
La misma información se puede consultar de manera más sintética.
st_dimensions(x)## from to offset delta refsys point values x/y
## x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [x]
## y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [y]
## band 1 6 NA NA NA NA NULLLectura y escritura de datos
Para leer y escribir datos desde archivos externos se utilizan las funciones read_stars y write_stars. Por ejemplo, carguemos un archivo que viene en el paquete stars:
filename <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(filename)
x## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## L7_ETMs.tif 1 54 69 68.91242 86 255
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [x]
## y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [y]
## band 1 6 NA NA NA NA NULLPara escribir un vector al disco:
write_stars(x,
"x_exp.tif")Las funciones st_read y st_write contienen más argumentos para definir ## Visualización de datos
Por último, el visualizar los datos nos puede dar una muy buena idea de los productos intermedios en un flujo de trabajo o verificar que la información que importamos o exportamos es la correcta. Para ver el objecto stars con todos sus atributos:
plot(x)
También se puede hacer un RGB en color natural o falso color.
plot(x,
rgb = c(3,2,1))
plot(x,
rgb = c(4,3,2))
Descargar información de trabajo
Primero vamos a cargar información desde mi github.
## [1] "D:\\Drive\\Jonathan_trabaggio\\Doctorado\\R\\SIG\\Curso\\Data\\Sentinel2-2A_10B_2020-01-01_2020-12-30.tif"
## [2] "D:\\Drive\\Jonathan_trabaggio\\Doctorado\\R\\SIG\\Curso\\Data\\SRTM_area.tif"
## [3] "D:\\Drive\\Jonathan_trabaggio\\Doctorado\\R\\SIG\\Curso\\Data\\roi.zip"Veamos qué descargamos.
im1 <- read_stars(paste0(getwd(), "/Data/", "Sentinel2-2A_10B_2020-01-01_2020-12-30.tif"))
im1## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2-2A_10B_2020-01-01_20... 6 913 1552 1501.178 2010.75 5041
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 103 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 81 19.6582 -0.000179663 WGS 84 FALSE NULL [y]
## band 1 10 NA NA NA NA B2,...,B12plot(im1)
Operaciones básicas
Elegir parte basado en los índices de los objetos stars
Para elegir parte de los rasters, se puede hacerlo eligiendo parte del arreglo. Esto se puede hacer utilizando la notación tradicional de objetos stars o usando la forma tidyverse-esque.
En la primera forma se utiliza la notación de arreglos [], en la que la primera dimensión corresponde a la dimensión de atributos, la segunda a la primera dimensión espacial (horizontal, x), la tercera a la segunda dimensión temporal (vertical, y) y la cuarta corresponde a la dimensión espectral (bandas).
im1[,1:10,1:20]## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2-2A_10B_2020-01-01_20... 73 691.25 1149 1315.034 1693.25 5041
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 10 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 20 19.6582 -0.000179663 WGS 84 FALSE NULL [y]
## band 1 10 NA NA NA NA B2,...,B12im1_sub1 <- im1[,1:10,1:20, 2]
im1_sub1## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2-2A_10B_2020-01-01_20... 189 485.75 556.5 555.88 635 889
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 10 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 20 19.6582 -0.000179663 WGS 84 FALSE NULL [y]
## band 2 2 NA NA NA NA B3library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x readr::parse_date() masks curl::parse_date()im1_sub2 <- im1 |>
slice(x, 1:10) |>
slice(y, 1:20) |>
slice(band, 2)
im1_sub2## stars object with 2 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2-2A_10B_2020-01-01_20... 189 485.75 556.5 555.88 635 889
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 10 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 20 19.6582 -0.000179663 WGS 84 FALSE NULL [y]Elegir parte basado en los valores de los objetos stars
Además podríamos estar interesados en filtrar valores de los rasters en función de los valores de cierta banda. Para ello, podemos hacer la función filter.
library(cubelyr)
im1_sub3 <- im1 |>
filter(x < -101.2332, x > -101.235,
y > 19.6564, y < 19.6582,
band > 3 )
plot(im1_sub3)
Elegir parte basado en un vector
Para seleccionar parte de un raster basado en un vector se puede utilizar la notación tradicional de corchetes [] o la función st_crop.
roi <- st_read(paste0(getwd(), "/Data/", "roi.shp"))## Reading layer `roi' from data source
## `D:\Drive\Jonathan_trabaggio\Doctorado\R\SIG\Curso\Data\roi.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -101.2307 ymin: 19.64597 xmax: -101.2206 ymax: 19.65094
## Geodetic CRS: WGS 84plot(roi)
im1_crop <- im1[roi]
plot(im1_crop)
im1_crop2 <- st_crop(im1,
roi)
plot(im1_crop2)Extraer parte de los objetos stars
Para extraer los valores del objeto stars y convertirlos a arreglos se puede usar la función pull. Usando esta función se extrae el atributo indicado dentro de pull.
im_pull <- im1 |>
pull(1)
# Mostrar únicamente el encabezado y algunos datos de muestra
head(im_pull)## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 434 482 414 401 414 415 386 333 200 205 227 264 259 317
## [2,] 434 482 441 398 401 411 440 513 533 478 372 398 397 411
## [3,] 443 487 400 354 339 337 371 399 479 494 462 490 444 367
## [4,] 488 441 326 309 325 333 359 324 390 511 501 489 463 393
## [5,] 560 409 309 320 314 244 114 84 104 303 466 538 402 367
## [6,] 594 418 321 306 255 73 134 225 296 300 208 453 471 275
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 155 316 402 424 426 434 422 403 369 367 189 456
## [2,] 244 224 357 419 382 398 420 392 350 346 270 408
## [3,] 192 171 257 422 478 415 394 347 336 268 459 251
## [4,] 149 153 160 343 456 465 397 320 219 193 247 109
## [5,] 130 133 136 292 416 447 435 372 171 211 155 125
## [6,] 115 107 100 207 348 460 444 373 169 559 194 225
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 350 625 290 457 886 977 748 984 730 728 587 415
## [2,] 268 1015 814 1119 807 531 406 378 347 405 436 324
## [3,] 317 317 561 603 642 625 819 470 448 579 712 886
## [4,] 259 278 664 718 374 946 814 365 285 474 713 661
## [5,] 199 328 831 1035 447 995 588 414 330 734 1551 998
## [6,] 308 485 1242 799 493 636 294 890 1186 1449 902 803
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 527 548 722 632 975 955 1417 828 645 1251 639 815
## [2,] 631 1060 1049 717 1335 667 679 596 724 508 521 738
## [3,] 431 304 688 513 662 932 549 375 632 611 535 555
## [4,] 506 383 581 582 633 1319 342 611 626 748 1037 477
## [5,] 298 279 393 629 716 1144 907 898 860 786 923 870
## [6,] 590 1003 1012 796 1155 1320 750 1144 1341 450 595 1233
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 612 304 149 105 198 254 783 756 711 627 495 680
## [2,] 965 639 185 110 174 234 817 699 891 907 597 560
## [3,] 535 444 477 204 130 275 484 974 749 890 679 692
## [4,] 690 424 1140 392 102 174 210 596 733 927 937 770
## [5,] 1019 825 662 505 321 89 100 477 405 1350 1211 770
## [6,] 618 818 827 673 420 160 100 477 405 1350 1211 1016
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 271 334 467 1429 396 781 1410 847 718 570 775 564
## [2,] 507 636 564 637 1280 1241 937 985 1004 783 690 1146
## [3,] 909 832 610 965 1300 1042 818 985 1004 783 690 1146
## [4,] 1098 1088 1220 965 1300 1042 818 928 1225 809 882 1280
## [5,] 1098 1088 1220 1130 713 881 1396 647 468 549 1516 1060
## [6,] 1078 756 856 620 444 460 1285 798 792 784 771 613
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1169 1099 1250 1012 1258 1114 776
## [2,] 1169 1099 1250 1012 1119 1352 1038
## [3,] 936 1132 821 704 883 937 1143
## [4,] 885 766 896 1327 1115 1047 582
## [5,] 707 1021 1393 1542 964 1324 844
## [6,] 558 809 1572 1467 1384 922 872
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 662 695 612 616 600 651 617 623 441 440 454 485 501 591
## [2,] 662 695 640 587 599 611 634 754 866 779 655 636 656 708
## [3,] 685 707 603 507 489 507 552 614 666 721 665 697 668 692
## [4,] 729 661 505 458 463 493 540 488 574 719 689 680 708 722
## [5,] 803 619 470 475 471 369 206 202 250 477 682 742 636 749
## [6,] 889 619 469 473 404 189 291 449 508 509 442 670 774 666
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 354 616 622 618 598 606 601 575 524 509 388 728
## [2,] 576 532 631 622 542 570 591 545 511 529 520 686
## [3,] 509 497 567 638 657 579 556 519 531 465 695 471
## [4,] 482 475 435 590 675 656 583 505 381 389 441 279
## [5,] 494 489 495 580 625 635 617 539 347 448 299 368
## [6,] 491 478 451 561 631 633 641 540 371 804 439 458
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 546 765 525 739 1014 1196 1030 1196 997 994 832 638
## [2,] 451 1211 1060 1333 966 722 666 633 580 655 739 563
## [3,] 548 527 811 892 931 873 1113 699 738 845 957 1136
## [4,] 475 556 897 1014 573 1112 1086 651 518 678 919 902
## [5,] 460 586 1025 1267 611 1300 825 651 537 1022 1935 1227
## [6,] 646 769 1521 1064 811 927 591 1132 1436 1634 1180 1050
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 681 855 957 846 1254 1153 1607 1181 1036 1885 869 1012
## [2,] 872 1281 1247 841 1392 965 1011 895 1094 819 830 1100
## [3,] 608 456 929 759 964 1185 758 657 850 854 856 903
## [4,] 785 649 765 836 752 1514 586 819 825 943 1472 665
## [5,] 545 540 611 863 900 1365 1111 1150 1021 1189 1111 1070
## [6,] 783 1256 1188 923 1295 1526 1052 1385 1653 757 971 1402
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 959 538 309 253 432 465 975 1008 885 806 816 906
## [2,] 1212 1011 456 320 361 477 1037 898 1030 1163 869 737
## [3,] 914 779 779 418 308 526 662 1172 1008 1130 937 864
## [4,] 830 760 1406 588 325 439 449 874 956 1112 1034 1040
## [5,] 1331 1105 985 770 577 320 275 693 605 1551 1448 1040
## [6,] 937 1052 1121 1041 665 396 275 693 605 1551 1448 1204
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 451 540 728 1702 667 1057 1646 1067 985 826 1027 897
## [2,] 697 734 791 996 1432 1421 1058 1169 1172 960 949 1411
## [3,] 1119 1122 809 1194 1556 1278 1054 1169 1172 960 949 1411
## [4,] 1283 1263 1537 1194 1556 1278 1054 1108 1492 1007 1002 1552
## [5,] 1283 1263 1537 1283 904 1088 1728 850 613 731 1772 1327
## [6,] 1251 1082 1037 850 673 695 1369 959 967 915 881 744
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1321 1295 1544 1302 1529 1336 996
## [2,] 1321 1295 1544 1302 1327 1576 1240
## [3,] 1129 1317 1073 879 1000 1232 1397
## [4,] 1023 1010 1173 1584 1252 1209 798
## [5,] 914 1130 1583 1844 1136 1448 1140
## [6,] 821 942 1729 1615 1701 1130 1033
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 939 986 839 840 824 881 815 791 414 445 507 580 556 644
## [2,] 939 986 886 780 817 857 892 1007 1127 969 784 759 824 799
## [3,] 960 978 815 696 674 693 765 842 941 993 898 944 898 782
## [4,] 1028 929 683 636 652 671 741 685 792 987 970 954 981 801
## [5,] 1091 847 653 671 656 505 295 227 272 559 905 1021 843 737
## [6,] 1203 866 662 671 563 226 277 420 589 595 453 922 976 569
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 367 625 815 858 834 805 787 778 728 705 443 767
## [2,] 553 524 727 820 764 782 818 764 701 708 554 730
## [3,] 454 397 610 865 898 781 773 732 737 601 738 483
## [4,] 353 330 360 726 895 884 805 676 533 414 423 272
## [5,] 282 263 286 624 833 857 854 758 408 404 311 320
## [6,] 230 232 226 408 778 864 878 744 416 869 422 457
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 609 830 522 732 1120 1357 1204 1245 1161 1185 1240 841
## [2,] 528 1335 1189 1544 1136 817 773 636 645 895 916 629
## [3,] 586 567 868 1100 1091 990 1242 726 788 1283 1275 1272
## [4,] 490 570 1190 1682 688 1233 1172 681 515 751 1027 1187
## [5,] 408 543 1135 1539 721 1415 914 717 593 1200 2212 1484
## [6,] 643 790 1792 1255 916 1064 593 1194 1500 1775 1272 1082
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 876 1056 1120 993 1360 1319 1735 1296 1241 2006 1024 1126
## [2,] 993 1421 1470 1018 1670 1098 1176 988 1367 1045 1038 1325
## [3,] 642 500 1034 885 1149 1408 877 790 990 1044 1079 1037
## [4,] 1145 756 896 899 885 1620 611 1069 985 1139 1517 997
## [5,] 573 521 609 912 1087 1671 1265 1341 1180 1239 1373 1230
## [6,] 800 1234 1248 1130 1539 1803 1195 1439 1708 977 1147 1449
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 1012 564 291 224 385 429 1113 1044 1026 1023 881 981
## [2,] 1635 1317 471 236 334 469 1140 1078 1100 1474 1019 857
## [3,] 952 966 907 357 272 516 719 1287 1128 1213 991 1103
## [4,] 924 835 1564 538 220 345 412 994 1045 1182 1213 1209
## [5,] 1444 1176 1157 883 553 217 235 688 719 1683 1504 1209
## [6,] 1078 1210 1274 1200 773 345 235 688 719 1683 1504 1300
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 455 525 749 1680 617 1122 1724 1225 1098 935 1173 1113
## [2,] 726 779 853 958 1472 1538 1244 1396 1339 1076 1079 1728
## [3,] 1242 1224 954 1333 1876 1471 1211 1396 1339 1076 1079 1728
## [4,] 1396 1567 1734 1333 1876 1471 1211 1232 1606 1202 1167 1777
## [5,] 1396 1567 1734 1416 1257 1240 1864 883 657 994 1970 1409
## [6,] 1354 1159 1164 1059 894 710 1565 1160 1098 983 945 831
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1560 1543 1869 1469 1919 1555 1183
## [2,] 1560 1543 1869 1469 1555 1816 1433
## [3,] 1374 1456 1202 985 1234 1628 1580
## [4,] 1195 1118 1310 1637 1316 1255 870
## [5,] 1039 1245 1744 1939 1295 1668 1282
## [6,] 874 1022 2013 1786 1940 1452 1289
##
## , , 4
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1078 1072 972 991 1036 1066 1187 1172 970 906 864 815 1065 984
## [2,] 1078 1072 1017 926 940 982 1002 1219 1378 1251 1072 1032 1153 1201
## [3,] 1096 1111 961 809 768 777 859 959 1105 1220 1109 1049 1115 1211
## [4,] 1168 1063 815 733 744 788 827 816 876 1075 1061 1058 1116 1265
## [5,] 1271 1011 752 759 770 566 360 409 508 913 1153 1147 1177 1222
## [6,] 1345 1149 759 783 610 399 597 850 984 867 856 1062 1404 1213
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 878 1081 1065 962 936 926 922 856 852 820 823 946
## [2,] 1080 1010 1130 980 883 879 920 882 846 786 1001 1044
## [3,] 1023 986 1057 1068 1006 926 891 847 803 860 1032 809
## [4,] 1036 924 903 1013 1090 986 899 784 658 732 669 551
## [5,] 1012 906 916 1030 1082 970 970 869 650 756 668 685
## [6,] 957 939 897 979 1104 1019 991 862 872 972 841 897
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 915 1058 1065 1217 1259 1300 1440 1491 1433 1402 1426 1152
## [2,] 828 1222 1472 1894 1248 1188 1202 1008 936 1345 1358 945
## [3,] 906 925 1156 1431 1159 1343 1558 1148 1110 1402 1446 1331
## [4,] 865 876 1507 1816 1132 1399 1410 1025 918 1027 1348 1651
## [5,] 920 877 1381 1591 1047 1581 1120 1009 1036 1536 2175 1631
## [6,] 1049 1166 1968 1552 1234 1511 1127 1468 1685 1813 1604 1104
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1152 1271 1342 1303 1393 1476 1599 1530 1598 2419 1517 1444
## [2,] 1137 1479 1437 1310 1760 1493 1534 1576 1671 1412 1397 1776
## [3,] 1102 949 1173 1106 1460 1397 1234 1114 1246 1394 1486 1492
## [4,] 1455 1055 1103 1050 1259 1439 1068 1205 1262 1323 1694 1368
## [5,] 994 932 1059 1036 1474 1834 1361 1465 1416 1593 1789 1551
## [6,] 1114 1353 1431 1287 1635 1705 1379 1710 1736 1419 1362 1583
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 1382 1005 628 618 711 829 1245 1359 1236 1280 1342 1248
## [2,] 1997 1713 965 636 671 856 1266 1281 1225 1513 1423 994
## [3,] 1379 1365 1353 749 641 846 1048 1339 1430 1199 1243 1381
## [4,] 1167 1271 1540 894 600 775 842 1190 1363 1310 1470 1378
## [5,] 1636 1532 1492 1257 924 643 585 891 1016 1775 1637 1378
## [6,] 1455 1319 1593 1579 1111 733 585 891 1016 1775 1637 1303
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 826 840 1089 1202 1057 1402 1680 1461 1422 1179 1425 1638
## [2,] 960 1158 1129 1368 1813 1809 1366 1402 1409 1307 1331 1762
## [3,] 1417 1342 1145 1510 1882 1604 1530 1402 1409 1307 1331 1762
## [4,] 1489 1673 1701 1510 1882 1604 1530 1368 1437 1413 1609 1731
## [5,] 1489 1673 1701 1474 1472 1514 1852 1304 805 1354 1913 1534
## [6,] 1446 1457 1437 1396 1131 1169 1565 1226 1181 1248 1116 1077
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1559 1759 1940 1692 1996 1627 1583
## [2,] 1559 1759 1940 1692 1660 1795 1639
## [3,] 1665 1445 1305 1273 1377 1759 1708
## [4,] 1301 1330 1417 1644 1461 1435 1161
## [5,] 1209 1386 1822 2042 1571 1782 1578
## [6,] 1114 1063 2397 2174 1617 1621 1438
##
## , , 5
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1200 1128 1139 1219 1298 1410 1644 1632 1762 1640 1350 1337 1648 1754
## [2,] 1200 1128 1067 1021 1025 1038 1151 1524 1966 1979 1676 1638 1694 2047
## [3,] 1179 1156 1077 975 886 883 928 1074 1225 1353 1207 1161 1203 1954
## [4,] 1237 1165 910 755 758 906 967 863 937 1141 1160 1182 1397 2180
## [5,] 1313 1106 802 827 814 547 502 704 1175 1507 1528 1286 1575 2381
## [6,] 1624 1225 871 834 617 676 1375 1866 1766 1522 1538 1523 1959 2747
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 1557 2128 1529 1209 1099 1021 951 968 1004 879 1350 1877
## [2,] 2076 1990 1679 1201 957 1001 1035 964 933 1049 1684 1594
## [3,] 2337 2208 1901 1501 1145 1003 994 902 958 1141 1752 1326
## [4,] 2777 2672 2099 1621 1295 1097 1021 876 759 1455 1315 1380
## [5,] 2996 3032 2832 2120 1414 1156 1116 869 1138 1458 1207 1850
## [6,] 3236 3008 3117 2547 1676 1212 1124 1000 1478 1959 1540 1930
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1329 1547 1798 2064 1537 1834 1965 1912 2004 1986 1857 1515
## [2,] 1355 1660 1991 2051 1457 1586 1967 1879 1566 1855 1930 1472
## [3,] 1458 1342 1868 2146 1433 1555 2115 1976 1692 1736 1850 1757
## [4,] 1778 1770 2042 2094 1233 1822 1913 1844 1733 1731 1842 1938
## [5,] 1983 2024 1909 1768 1421 2012 1754 1626 1603 2099 2435 2011
## [6,] 2058 1856 2177 1445 1991 1967 1555 1910 1965 1904 1679 1465
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1462 1630 1485 1627 1912 1767 2125 2427 2285 2109 1716 2190
## [2,] 1486 1650 1721 1634 1959 1857 1946 2082 2135 2136 2199 2470
## [3,] 1690 1542 1368 1705 1893 1805 1850 1923 1693 1850 2163 2432
## [4,] 1928 1698 1434 1519 1644 2051 1533 1601 1569 1955 2314 1829
## [5,] 1718 1741 1547 1610 1638 1888 1591 1679 1968 2188 1980 1913
## [6,] 1562 1847 1877 1686 1854 1804 1520 1944 2213 2000 2193 2228
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 2382 1954 1558 1428 1670 1483 1504 1853 1591 1747 2130 1765
## [2,] 2567 2400 1671 1514 1624 1708 1754 1746 1498 1629 1896 1639
## [3,] 2348 2013 1803 1597 1469 1470 1536 1653 1748 1517 1812 1592
## [4,] 1810 2024 1959 1585 1874 1776 1464 1664 1708 1615 1507 1403
## [5,] 1982 2037 2288 1939 1946 1954 1510 1601 1493 1891 1844 1403
## [6,] 1932 1907 2137 2202 1881 1470 1510 1601 1493 1891 1844 1577
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1484 1494 1841 2080 1903 1949 1799 1545 1717 1698 1758 2023
## [2,] 1602 1593 1639 1717 2291 2085 1667 1496 1505 1652 1732 2013
## [3,] 1429 1390 1426 1691 2100 1855 1655 1496 1505 1652 1732 2013
## [4,] 1618 1814 1851 1691 2100 1855 1655 1546 1657 1779 1748 1888
## [5,] 1618 1814 1851 1720 1623 1717 1959 1261 1362 1445 1896 1837
## [6,] 1759 1991 1862 1754 1892 1854 1958 1362 1493 1507 1589 1496
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1791 1780 1967 1714 2059 1715 1565
## [2,] 1791 1780 1967 1714 1740 1913 1757
## [3,] 1672 1807 1620 1514 1510 1841 1891
## [4,] 1419 1402 1459 1802 1704 1527 1646
## [5,] 1411 1494 2038 2236 1472 1757 1779
## [6,] 1466 1433 2122 2400 1977 1723 1572
##
## , , 6
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1246 1271 1268 1252 1259 1473 1689 2021 2088 1716 1603 1578 1866 1810
## [2,] 1246 1271 1156 1135 1156 1225 1301 1650 2214 2454 1938 1629 1776 2510
## [3,] 1244 1235 1187 1030 944 978 1039 1153 1344 1303 1184 1299 1457 2308
## [4,] 1314 1255 998 859 884 918 873 747 1009 1297 1332 1333 1479 2470
## [5,] 1423 1241 887 901 892 707 712 1070 1440 1638 1641 1412 1781 2882
## [6,] 1863 1350 918 883 700 700 1261 1945 1915 1928 1918 1800 2230 3187
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 1849 2326 1596 1186 1118 1198 1183 1094 985 886 1602 2139
## [2,] 2513 2332 1912 1361 1074 1070 1115 1085 1076 1230 1906 1696
## [3,] 2616 2477 2191 1648 1259 1132 1090 1010 1003 1350 2138 1628
## [4,] 3260 3198 2465 1807 1460 1258 1149 972 951 1509 1487 1491
## [5,] 3622 3610 3435 2520 1605 1267 1192 1019 1316 1717 1537 2025
## [6,] 3726 3651 3742 3216 2019 1421 1251 1125 1738 1966 1727 2234
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1511 1755 2200 2470 1555 1677 2148 2163 2061 2064 1985 1577
## [2,] 1314 1493 2081 2229 1594 1898 2348 2029 1793 2051 2071 1670
## [3,] 1902 1697 2160 2181 1519 1735 2427 2284 1865 1762 1962 1766
## [4,] 1982 1933 2153 2228 1323 1829 2008 2251 1918 1857 1944 2165
## [5,] 2367 2246 2015 1732 1535 2185 1723 1673 1859 2271 2435 2006
## [6,] 2282 2113 2529 1546 2293 2156 1850 2072 2081 2026 1863 1567
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1680 1821 1698 1757 1837 1798 2103 2549 2459 2456 2008 2468
## [2,] 1785 1704 1714 1712 2101 1942 2213 2438 2406 2421 2514 2760
## [3,] 1641 1601 1542 1860 2115 1928 2163 2077 1771 2085 2452 2685
## [4,] 2326 2009 1676 1669 1710 1909 1681 1809 1758 2125 2498 2038
## [5,] 1778 1892 1823 1670 1866 2107 1620 1760 2172 2375 2149 2013
## [6,] 1795 1866 1881 1803 1866 1742 1635 2080 2371 2287 2324 2391
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 2667 2159 1796 1857 1926 1716 1678 2097 1639 1949 2383 1980
## [2,] 2757 2579 1936 1843 1842 1841 1832 1826 1487 1778 2028 1818
## [3,] 2673 2321 2038 1765 1694 1765 1747 1830 1838 1593 1901 1717
## [4,] 2032 2332 2072 1836 2172 1978 1713 1849 1816 1669 1576 1407
## [5,] 2112 2324 2539 2286 2254 2290 1777 1797 1634 2097 1866 1407
## [6,] 2103 1920 2235 2403 2127 1853 1777 1797 1634 2097 1866 1614
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1673 1724 2010 2182 2190 2086 1861 1566 1898 1832 1918 2167
## [2,] 1786 1784 1879 1911 2436 2204 1747 1454 1632 1794 1878 2032
## [3,] 1563 1437 1525 1782 2212 1886 1744 1454 1632 1794 1878 2032
## [4,] 1636 1829 1840 1782 2212 1886 1744 1590 1718 1803 1827 1949
## [5,] 1636 1829 1840 1767 1674 1811 1994 1409 1470 1512 1957 1918
## [6,] 1894 2131 1937 1895 2030 2091 1997 1426 1709 1636 1725 1634
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1768 1820 1934 1744 2073 1771 1707
## [2,] 1768 1820 1934 1744 1716 1966 1873
## [3,] 1749 1821 1744 1660 1576 1896 1851
## [4,] 1480 1518 1515 1967 1706 1693 1798
## [5,] 1452 1554 2137 2192 1496 1802 1790
## [6,] 1653 1398 2207 2378 2023 1814 1656
##
## , , 7
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1318 1372 1324 1338 1383 1619 1868 2271 2197 1884 1620 1586 1925 2182
## [2,] 1318 1372 1248 1233 1357 1336 1343 1686 2531 2585 2089 1937 2015 2760
## [3,] 1327 1354 1267 1120 1005 1046 1134 1208 1398 1474 1348 1417 1455 2533
## [4,] 1386 1399 1104 938 976 989 1095 1012 1218 1463 1427 1424 1619 2743
## [5,] 1504 1322 973 1010 960 807 591 912 1249 1751 1739 1581 1872 3018
## [6,] 2108 1390 1028 986 858 668 1523 2383 2210 1974 2039 1880 2463 3438
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 1877 2681 1633 1361 1305 1241 1273 1231 1191 1184 1522 2377
## [2,] 2719 2481 2175 1432 1175 1184 1238 1179 1156 1303 2172 1665
## [3,] 2896 2767 2424 1811 1404 1219 1186 1135 1177 1373 2353 1627
## [4,] 3585 3357 2545 2059 1570 1396 1225 1061 950 1710 1621 1626
## [5,] 3939 3867 3723 2707 1707 1370 1350 1208 1208 2090 1312 2346
## [6,] 4087 3874 3986 3524 2174 1500 1415 1348 1567 2469 1816 2279
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1722 1797 2202 2239 1625 1780 2396 2099 2344 2174 1848 1830
## [2,] 1394 1892 1849 2119 1670 1984 2374 2222 1653 2178 2127 1737
## [3,] 2032 1381 2425 2169 1562 1873 2507 2394 2218 1848 2035 1874
## [4,] 2115 2123 2305 2324 1226 1993 2017 2457 1888 1952 2005 2350
## [5,] 2491 2420 2121 1790 1553 2287 1843 1886 1883 2311 2574 2068
## [6,] 2684 2239 2320 1481 2620 1956 1845 2195 2018 2092 1816 1630
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1533 1824 1693 1848 1899 1742 2413 2753 2602 2981 1949 2619
## [2,] 1928 1861 1846 1736 2265 1940 2269 2425 2705 2716 2677 2977
## [3,] 1727 1735 1588 1955 2163 1986 2291 2329 1895 2345 2800 2914
## [4,] 2383 2173 1646 1861 1613 2012 1922 1935 1800 2097 2588 2047
## [5,] 1881 2161 1770 1796 1900 2081 1680 1802 2199 2640 2311 2141
## [6,] 1799 1941 2070 1960 1962 1975 1612 2234 2312 2300 2366 2544
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 2912 2447 1906 1811 2109 1852 1725 2124 1716 2052 2439 2086
## [2,] 2773 2819 1905 2007 1927 2062 1967 2083 1575 1883 2370 1908
## [3,] 2931 2385 2157 2057 1724 2038 1586 1912 1906 1542 1821 1644
## [4,] 2016 2653 2424 1837 2543 2179 1829 1932 1859 1625 1687 1584
## [5,] 2144 2337 2860 2168 2335 2223 1657 1984 1669 2251 2015 1584
## [6,] 2292 1955 2405 2697 2191 2222 1657 1984 1669 2251 2015 1617
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1795 1895 2036 2427 2226 2203 2009 1459 2107 1832 2068 2095
## [2,] 1908 1801 1875 1840 2575 2264 1720 1454 1803 1848 1940 2218
## [3,] 1629 1686 1621 1881 2240 1870 1685 1454 1803 1848 1940 2218
## [4,] 1655 1770 1962 1881 2240 1870 1685 1472 1792 1778 1754 2105
## [5,] 1655 1770 1962 1792 1788 1796 2444 1106 1696 1553 2261 1937
## [6,] 2043 2077 1919 1936 2194 2040 2060 1429 1828 1675 1594 1778
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1719 1926 2236 1721 2040 1912 1753
## [2,] 1719 1926 2236 1721 1778 1964 1715
## [3,] 1568 1858 1754 1712 1537 1862 1890
## [4,] 1493 1620 1713 2018 1583 1784 1802
## [5,] 1422 1692 1954 2276 1539 1938 1585
## [6,] 1583 1521 2389 2206 2269 1672 1637
##
## , , 8
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1412 1394 1445 1371 1454 1671 1929 2432 2313 1977 1834 1773 2173 2053
## [2,] 1412 1394 1294 1263 1382 1416 1473 1878 2560 2673 2160 1811 2094 2717
## [3,] 1372 1435 1330 1132 1076 1088 1196 1285 1531 1485 1428 1431 1725 2587
## [4,] 1380 1393 1172 1022 1052 1064 993 957 1185 1500 1455 1454 1691 2865
## [5,] 1588 1406 1041 1029 1031 808 756 1198 1447 1799 1861 1607 2056 3175
## [6,] 2062 1584 1034 1059 839 727 1419 2113 2044 2031 2054 1991 2653 3539
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 2173 2465 1946 1403 1307 1403 1407 1279 1145 1155 1779 2189
## [2,] 2870 2581 2235 1543 1270 1246 1268 1283 1285 1316 2134 1908
## [3,] 2947 2770 2475 1856 1500 1310 1225 1153 1158 1611 2249 1724
## [4,] 3585 3511 2702 2131 1731 1451 1351 1119 1079 1680 1703 1672
## [5,] 3940 4064 3669 2818 1843 1493 1433 1248 1419 2061 1688 2169
## [6,] 4042 4030 3999 3535 2313 1629 1429 1387 1898 2101 2118 2352
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1746 1752 2485 2700 1604 1841 2367 2204 2198 2242 1964 1722
## [2,] 1366 1642 2174 2475 1711 2114 2605 2280 1961 2263 2122 1891
## [3,] 2057 1842 2349 2311 1626 1899 2563 2468 1898 1886 2040 1827
## [4,] 2106 1976 2351 2207 1484 1797 2201 2406 2085 2083 2084 2380
## [5,] 2740 2355 2149 1681 1683 2279 1860 1968 1978 2392 2564 2084
## [6,] 2477 2429 2368 1963 2235 2166 1964 2153 2092 2094 1939 1630
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1769 1825 1868 1863 1897 1956 2346 2625 2677 2665 2085 2741
## [2,] 1918 1813 1839 1874 2228 2031 2388 2593 2850 2725 2817 3007
## [3,] 1834 1772 1680 2028 2229 1964 2407 2426 2047 2384 2821 2931
## [4,] 2506 2224 1858 1856 1779 1823 1963 1960 1900 2286 2628 2162
## [5,] 1957 2142 2044 1810 1999 2038 1774 1861 2187 2505 2284 2189
## [6,] 1801 1956 2050 1895 1948 1935 1770 2205 2367 2481 2386 2605
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 2926 2511 2102 1960 2096 1942 1772 2162 1758 2006 2476 2195
## [2,] 2950 2677 2325 2102 1964 2131 1977 2079 1559 1973 2229 2011
## [3,] 2984 2604 2235 2055 1812 1996 1744 1980 1857 1740 1984 1798
## [4,] 2116 2613 2279 1974 2260 2268 1807 1993 1916 1713 1691 1515
## [5,] 2146 2526 2706 2472 2441 2412 1955 1925 1798 2126 1975 1515
## [6,] 2291 2102 2500 2582 2298 2064 1955 1925 1798 2126 1975 1714
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1831 1911 2060 2307 2308 2268 1876 1717 2042 1854 2076 2137
## [2,] 1860 1934 1887 1902 2428 2438 1733 1562 1807 1868 2019 2064
## [3,] 1630 1571 1666 1891 2216 1942 1798 1562 1807 1868 2019 2064
## [4,] 1749 1849 1920 1891 2216 1942 1798 1539 1722 1788 1879 2041
## [5,] 1749 1849 1920 1824 1707 1929 2184 1502 1370 1688 2069 1960
## [6,] 2026 2112 1967 2005 2113 2306 2022 1699 1742 1769 1813 1854
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 1694 1868 2024 1818 2018 1864 1946
## [2,] 1694 1868 2024 1818 1790 1931 1953
## [3,] 1837 1783 1862 1739 1682 1746 1894
## [4,] 1620 1607 1680 1926 1647 1846 1854
## [5,] 1548 1698 2210 2162 1637 1869 1679
## [6,] 1787 1533 2376 2222 2003 1912 1755
##
## , , 9
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1909 1893 1870 1809 1834 1967 1998 2006 1730 1550 1696 1754 1978 1788
## [2,] 1909 1893 1742 1742 1818 1865 1990 2189 2252 2194 1940 1885 2054 2223
## [3,] 1856 1882 1761 1668 1603 1605 1675 1833 2036 2094 2033 2083 2098 2244
## [4,] 1908 1849 1679 1559 1533 1538 1469 1404 1612 2009 2050 2091 2134 2229
## [5,] 2054 1851 1587 1568 1497 1185 928 913 1108 1532 1976 2131 2321 2226
## [6,] 2588 2151 1665 1508 1219 900 1004 1328 1568 1577 1647 2057 2504 2082
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 1584 2019 2279 2168 2086 2095 2080 2036 1950 1803 1687 1695
## [2,] 2037 2000 2217 2160 1997 1959 1956 1895 1844 1777 1760 1584
## [3,] 1911 1791 2130 2346 2249 2050 1929 1831 1780 1747 1726 1415
## [4,] 1875 1692 1831 2186 2433 2312 2056 1735 1493 1414 1303 1273
## [5,] 1777 1642 1701 2111 2347 2329 2189 1908 1516 1371 1220 1353
## [6,] 1596 1594 1631 1942 2329 2379 2279 2151 1900 1720 1534 1655
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1581 1630 1687 2141 2019 1978 2153 2314 2345 2077 1961 1870
## [2,] 1513 2000 2112 2465 1958 1924 2103 2029 1992 2082 1964 1726
## [3,] 1686 1609 2043 2398 1857 1863 2132 1926 1939 2156 2089 1906
## [4,] 1679 1734 2225 2376 1765 1947 2030 1862 1637 1916 2276 2344
## [5,] 1708 1807 2121 2133 1764 2301 1733 1573 1747 2408 2886 2476
## [6,] 1911 1980 2384 1798 1876 1955 1642 1790 2202 2453 2443 2018
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1976 2117 2086 2009 2063 1994 2353 2468 2635 2651 2375 2556
## [2,] 1835 2139 2261 2004 2174 1985 2145 2540 3014 2929 2730 2865
## [3,] 1664 1697 2153 2061 2120 2165 2145 2427 2581 2638 2795 2660
## [4,] 2209 1961 2043 1854 1840 2049 1936 2202 2110 2286 2362 2143
## [5,] 1935 1682 1786 1807 1918 2176 2159 2269 2195 2266 2239 2193
## [6,] 1775 1966 2067 2012 2381 2448 2216 2347 2430 2430 2495 2429
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 2587 2065 1332 1106 1370 1454 1830 2025 1959 2113 2266 2046
## [2,] 2983 2730 1735 1174 1298 1531 2064 2154 1958 2042 2153 1914
## [3,] 2498 2520 2133 1381 1239 1494 1808 2077 2132 2050 2175 2010
## [4,] 2135 2292 2308 1610 1270 1386 1446 1778 1989 1998 1952 1990
## [5,] 2239 2307 2330 2017 1618 1232 1135 1481 1704 2258 2182 1990
## [6,] 2309 2227 2296 2356 1934 1372 1135 1481 1704 2258 2182 2049
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1600 1594 1784 1854 1849 2143 2227 2079 2055 2006 2024 2189
## [2,] 1676 1703 1668 1775 2124 2265 2141 2043 2022 2003 2116 2352
## [3,] 1947 1918 1890 2031 2403 2241 2077 2043 2022 2003 2116 2352
## [4,] 2202 2284 2188 2031 2403 2241 2077 1887 1886 1966 2103 2250
## [5,] 2202 2284 2188 2085 2112 2250 2214 1558 1410 1807 2126 2086
## [6,] 2145 2119 2019 2012 2006 2052 2184 1826 1878 1935 1939 1875
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 2192 2409 2448 2288 2507 2349 2311
## [2,] 2192 2409 2448 2288 2330 2589 2409
## [3,] 2303 2215 2041 1864 2130 2585 2343
## [4,] 1970 1995 1911 2019 2000 2085 2116
## [5,] 1952 2167 2746 2542 2156 2304 2167
## [6,] 1868 1990 2826 2773 2373 2208 2178
##
## , , 10
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1480 1488 1448 1391 1399 1439 1500 1418 1118 1129 1194 1196 1329 1207
## [2,] 1480 1488 1411 1362 1404 1457 1524 1705 1699 1508 1438 1516 1524 1417
## [3,] 1469 1495 1384 1305 1276 1294 1363 1502 1696 1853 1784 1709 1639 1443
## [4,] 1483 1444 1331 1270 1258 1278 1250 1186 1353 1675 1802 1771 1624 1419
## [5,] 1602 1485 1308 1293 1249 1011 698 600 720 1149 1613 1787 1682 1344
## [6,] 1925 1700 1348 1263 1025 662 669 857 1100 1082 1161 1473 1732 1237
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 967 1219 1648 1758 1693 1707 1708 1682 1574 1444 1168 1107
## [2,] 1176 1180 1553 1760 1696 1647 1687 1609 1517 1366 1247 1223
## [3,] 1082 1017 1325 1785 1925 1795 1625 1540 1442 1296 1233 989
## [4,] 1006 831 1022 1517 1973 2025 1771 1460 1152 968 842 737
## [5,] 906 780 888 1302 1806 1975 1869 1597 1061 903 796 796
## [6,] 783 735 765 1068 1622 1921 1914 1711 1286 1166 1038 1002
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1197 1371 1348 1813 1784 1554 1666 1832 1869 1615 1513 1524
## [2,] 1164 1694 1944 2578 1904 1628 1507 1426 1576 1604 1460 1288
## [3,] 1116 1266 1580 1960 1555 1542 1636 1344 1440 1709 1716 1720
## [4,] 1070 1179 1772 2129 1641 1746 1638 1317 1138 1466 1891 1954
## [5,] 1035 1191 1759 2042 1479 1958 1379 1139 1276 2058 2795 2263
## [6,] 1151 1411 2318 2014 1456 1641 1171 1317 1878 2246 2229 1803
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1676 1897 1839 1746 1938 1834 1997 1900 1911 2084 1844 1859
## [2,] 1403 1894 2222 1769 1952 1768 1752 1898 2240 2094 1965 2159
## [3,] 1473 1427 1895 1624 1782 1838 1541 1656 1951 1920 1826 1863
## [4,] 1631 1315 1726 1651 1612 1829 1597 1779 1719 1817 1846 1722
## [5,] 1514 1261 1442 1531 1644 1984 2042 1953 1814 1855 1826 1836
## [6,] 1552 1746 1849 1825 2226 2318 2029 2183 2255 1993 1913 1926
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 1801 1435 821 613 811 1012 1657 1674 1638 1808 1789 1556
## [2,] 2256 1933 1138 682 757 1033 1753 1802 1710 1783 1756 1467
## [3,] 1830 1766 1574 936 724 971 1456 1813 1934 1876 1982 1869
## [4,] 1696 1680 1837 1220 727 833 949 1362 1823 1771 1810 1848
## [5,] 2060 1874 1744 1530 1115 746 692 1057 1369 1878 1919 1848
## [6,] 1903 1872 1928 1790 1474 860 692 1057 1369 1878 1919 1892
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1109 1173 1383 1471 1400 1760 2154 2059 1864 1707 1748 1814
## [2,] 1322 1417 1332 1466 1815 2213 2112 2035 1936 1818 1871 2119
## [3,] 1822 1868 1696 1806 2438 2347 2042 2035 1936 1818 1871 2119
## [4,] 2092 2214 2123 1806 2438 2347 2042 1855 1795 1812 1876 2060
## [5,] 2092 2214 2123 1972 1865 1921 1924 1386 1206 1498 1880 1915
## [6,] 2027 1868 1807 1722 1528 1570 1869 1689 1625 1704 1720 1563
## [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## [1,] 2026 2274 2359 2254 2495 2443 2402
## [2,] 2026 2274 2359 2254 2278 2625 2568
## [3,] 2199 2084 1830 1685 1920 2164 2106
## [4,] 1878 1942 1804 1891 1908 2030 1852
## [5,] 1824 2048 2751 2488 2102 2182 2015
## [6,] 1647 1812 2697 2733 2209 2019 2068Agregar un nuevo atributo
Para agregar un nuevo atributo en función de los que ya existen.
names(im1) <- "Sentinel2"
im2 <- im1 |>
mutate(band2 = 2 * Sentinel2)
im2## stars object with 3 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2 6 913 1552 1501.178 2010.75 5041
## band2 12 1826 3104 3002.356 4021.50 10082
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 103 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 81 19.6582 -0.000179663 WGS 84 FALSE NULL [y]
## band 1 10 NA NA NA NA B2,...,B12Ahora ya tenemos dos atributos.
Para seleccionar ciertos atributos
Para seleccionar algunos atributos se usa select.
im2 |>
select(band2)## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## band2 12 1826 3104 3002.356 4021.5 10082
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 103 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 81 19.6582 -0.000179663 WGS 84 FALSE NULL [y]
## band 1 10 NA NA NA NA B2,...,B12Cálculo de índices
En algunos casos estaremos más interesados en calcular índices que resalten algunas características de la superficie terrestre. Por ejemplo, algún índice de vegetación. Esto se puede realizar al seleccionar las bandas de interés y aplicar la fórmula.
NIR <- im1[,,,8]
R <- im1 [,,,4]
# Tidyverse-esque
# Si no se agrega drop = F, se elimina la tercera dimensión y luego ya no permite concatenar esta banda con el resto de la imagen.
NIR <- im1 |>
slice(band, 8, drop = F)
R <- im1 |>
slice(band, 4, drop = F)
NDVI <- (NIR - R) / (NIR + R)
plot(NDVI)
Agregar bandas a un raster
Para agregar una banda a un raster se utiliza la función c.
im1_cNDVI <- c(im1,
NDVI,
along = 3)Transformaciones de formato
Transformación de datos stars a sf
Para transformar un vector a raster se utiliza la función st_rasterize. Para ello, hay que determinar la resolución a la que se quiere rasterizar esta información. Utilizamos los valores de resolución que tiene el objeto im1. Recordemos que estos valores están en las unidades del CRS de la capa. En este caso, el CRS está en coordenadas geográficas, por eso se utilizarán valores en grados.
roi_rast <- st_rasterize(roi["id"], dx = 0.000179663, dy = -0.000179663)
plot(roi_rast)## Warning in plot.stars(roi_rast): breaks="quantile" leads to a single class;
## maybe try breaks="equal" instead?
Transformación de raster a vector
Para transformar de raster a un vector de puntos se puede utilizar la función st_as_sf, utilizando el argumento de as_points. Además, si se desean obtener los datos en formato largo en lugar de ancho se puede utilizar el argumento long.
st_as_sf(im1,
as_points = TRUE,
merge = FALSE)## Simple feature collection with 8343 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -101.2352 ymin: 19.64373 xmax: -101.2169 ymax: 19.6581
## Geodetic CRS: WGS 84
## First 10 features:
## B2 B3 B4 B5 B6 B7 B8 B8A B11 B12 geometry
## 1 434 662 939 1078 1200 1246 1318 1412 1909 1480 POINT (-101.2352 19.6581)
## 2 434 662 939 1078 1200 1246 1318 1412 1909 1480 POINT (-101.235 19.6581)
## 3 443 685 960 1096 1179 1244 1327 1372 1856 1469 POINT (-101.2348 19.6581)
## 4 488 729 1028 1168 1237 1314 1386 1380 1908 1483 POINT (-101.2347 19.6581)
## 5 560 803 1091 1271 1313 1423 1504 1588 2054 1602 POINT (-101.2345 19.6581)
## 6 594 889 1203 1345 1624 1863 2108 2062 2588 1925 POINT (-101.2343 19.6581)
## 7 461 693 980 1257 1480 1604 1709 1862 2636 1926 POINT (-101.2341 19.6581)
## 8 420 645 907 1067 1152 1321 1561 1701 2432 1814 POINT (-101.2339 19.6581)
## 9 436 635 920 1047 1187 1349 1544 1679 2429 1796 POINT (-101.2338 19.6581)
## 10 345 528 745 990 1185 1381 1333 1621 2329 1695 POINT (-101.2336 19.6581)Utilizando la misma función se puede convertir de raster a vector de polígonos. En este caso, se indica como argumentos que no se desea obtener una capa de puntos y en merge se puede indicar si se desea fusionar los polígonos con un mismo valor o no.
im1_poly <- st_as_sf(im1,
as_points = F,
merge = T)
plot(im1_poly)
Operaciones espaciales
Reproyección de rasters
Para reproyectar un raster se utiliza la función st_transform. Al igual que en sf se necesita indicar el CRS objetivo ya sea mediante el código EPSG o obteniendo dicha clave con st_crs. Este tipo de reproyección no presenta pérdida de información debido a las diferencias entre proyecciones. Por lo tanto, los pixeles pueden no ser de dimensiones homogeneas en toda la imagen. Ve la siguiente función para hacer reproyecciones más clásicas.
im1_utm <- st_transform(im1,
32614)
im1_utm## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Sentinel2 6 913 1552 1501.178 2010.75 5041
## dimension(s):
## from to offset delta refsys point
## x 1 103 NA NA WGS 84 / UTM zone 14N FALSE
## y 1 81 NA NA WGS 84 / UTM zone 14N FALSE
## band 1 10 NA NA NA NA
## values x/y
## x [103x81] 265619,...,267562 [x]
## y [103x81] 2173570,...,2175186 [y]
## band B2,...,B12
## curvilinear gridResampling
Para cambiar la resolución de un raster o reproyectar con pérdida de información se utiliza la función st_warp. Esta función permite definir el tamaño de celda objetivo mediante el argumento cellsize, así como el CRS objetivo. Esta función hace la reproyección donde todos los tamaños de celda son iguales y quizás es con la que cualquier usuario esté más familiarizado.
im1_utm_20m <- st_warp(im1,
crs = st_crs(im1),
cellsize = c(20, 20))
im1_utm_20m## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## Sentinel2 NA NA NA NaN NA NA 10
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 1 -101.235 20 WGS 84 NA NULL [x]
## y 1 1 19.6582 -20 WGS 84 NA NULL [y]
## band 1 10 NA NA NA NA B2,...,B12Crear mosaicos a partir de varios rasters
Para crear un mosaico a partir de dos o más imágenes se utiliza la función st_mosaic.
im_mosaic <- st_mosaic(im1_sub2, im1_sub3)
plot(im_mosaic)
Operaciones basadas en los valores de raster
Reclasificar
Para reclasificar un raster se puede realizar con la función cut, indicando los intervalos para las clases.
reclass_NDVI <- cut(NDVI,
breaks = c(-1, 0, 0.2, 0.5, 1),
labels = F)
plot(reclass_NDVI)
Enmascarar
Para enmascarar algún raster utilizando otro. De igual manera se pueden emascarar los valores que no cumplan con algún criterio de interés.
im_mask <- reclass_NDVI
im_mask[im_mask != 2] <- NA
reclass_NDVI2 <- reclass_NDVI
reclass_NDVI2[is.na(im_mask)] <- NA
plot(reclass_NDVI2)## Warning in plot.stars(reclass_NDVI2): breaks="quantile" leads to a single class;
## maybe try breaks="equal" instead?
Cálculo de métricas por banda
Para calcular métrias por banda se puede utilizar la función pull para obtener el atributo del que se desea obtener la información. En este caso la imagen solo cuenta con un atributo, así que seleccionamos esa entrada. Después, podemos utilizar la función apply para calcular la media sobre la dimensión 3 del arreglo. Otra forma de acceder a los datos en formato de arreglo es con [[]].
# Opción con pull
im1 |>
pull(1) |>
apply(3, mean)## [1] 527.6186 772.9070 935.5342 1230.8505 1700.1954 1874.5490 2004.2371
## [8] 2076.1721 2236.3476 1653.3689# Opción con doble corchete
im1[[1]] |>
apply(3, mean)## [1] 527.6186 772.9070 935.5342 1230.8505 1700.1954 1874.5490 2004.2371
## [8] 2076.1721 2236.3476 1653.3689Aplicar una función sobre todos los pixeles
Si se desea aplicar una función sobre todas las bandas de una imagen por pixel, se puede utilizar la función st_apply para facilitar el cálculo de variables. Esto es muy utilizado para análisis con series de tiempo y obtener por ejemplo la media de cada pixel.
mean_allBands <- st_apply(im1,
1:2,
mean)
mean_allBands## stars object with 2 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## mean 225.8 1253 1502.2 1501.178 1716.45 4223.5
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 103 -101.235 0.000179663 WGS 84 FALSE NULL [x]
## y 1 81 19.6582 -0.000179663 WGS 84 FALSE NULL [y]plot(mean_allBands)
Tabla de frecuencias
Para obtener una tabla de las frecuencias por el valor de pixel se utiliza la función table.
table(reclass_NDVI)## reclass_NDVI
## 1 2 3 4
## 64 2364 5527 388Objetos stars proxy
En algunos casos los archivos raster con los que se desea trabajar son muy pesados y no caben en la memoria disponible RAM. Por ello, al llamar read_stars se crea automáticamente un objeto proxy stars. Estos objetos permite cargar los datos únicamente cuando se vayan a utilizar. Para cargar los datos ya que se vayan a utilizar se usa la función st_as_stars.
Visualizacion avanzada
Para generar visualizaciones más avanzadas podemos utilizar otros paquetes como ggplot2 o tmap. Aquí veremos un ejemplo de cada uno.
ggplot2
En ggplot se puede determinar directamente el color y relleno de cada capa vectorial. En este ejemplo se carga la biblioteca viridis para usar dicha escala de color.
library(viridis)## Loading required package: viridisLite## Loading required package: viridisLite
ggplot() +
geom_stars(data = im1) +
coord_equal() +
facet_wrap(~band) +
theme_void() +
scale_fill_viridis() +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0))
tmap
La opción de tmap require de utilizar algún atributo de la información para determinar el color de relleno de cada polígono de acuerdo a los valores de ese atributo. Es similar a ggplot utilizando la opción de aes. Podemos ver algunas de las paletas que ya vienen pre hechas tanto para ggplot como para tmap con tmaptools::palette_explorer()
library(tmap)
# Cargar shapr
tm_shape(im1) +
# Elegir forma de visualización
tm_raster()
</div>
</div>