Geração de binários geogrid de uso da terra do MODIS para o WRF
Jonatan Tatsch
August 25, 2016
Delimitação do estado do RS
states <- readRDS("../../atividade1/data/br_states.rds")
pol_rs <- states[states@data$NAME_1 == "Rio Grande do Sul",]Cobertura da terra do MODIS
lct <- raster("../../atividade1/data/LCType.tif")
lct## class : RasterLayer
## dimensions : 43200, 86400, 3732480000 (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /home/hidrometeorologista/UFSM/ensino/_2semestre2016/isa/norton/atividade1/data/LCType.tif
## names : LCType
## values : 0, 255 (min, max)ll <- projection(lct)
ll## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"Tiles de 10 x 10 graus para o globo
r <- raster(ncols = 1200, nrows = 1200,
xmn = -60, xmx = -50,
ymx = -20, ymn = -30)
r## class : RasterLayer
## dimensions : 1200, 1200, 1440000 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -60, -50, -30, -20 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0Tiles que abrangem o RS
# raster de 10 x 10 graus
ref <- raster(res = 10)
ref[] <- 1:ncell(ref)
plot(ref)
plot(ref)
cels <- unlist(extract(ref, pol_rs))
cels## [1] 445 410 409# supondo que fosse só um tile
#cel1 <- cels[1]
# rcels <- ref %in% cel1
rcels <- ref %in% cels
plot(rcels, xlim = c(-180, 180))
plot(pol_rs, add = TRUE, col = "red")
grid(nx=36,ny=18)
text(ref, cex = 0.5)
Shape FEPAM
pol <- shapefile("../../atividade1/data/areas_umidas_FZB_LLsirgas2000.shp")
pol## class : SpatialPolygonsDataFrame
## features : 1703
## extent : -57.64514, -49.72959, -33.69034, -28.08131 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## variables : 2
## names : Hectares, NOME
## min values : 1.000082e+02, Banhado da Foz do Arroio Pastoreio
## max values : 9.991974e+00, Sistema Lagoa Pequenaprojection(pol)## [1] "+proj=longlat +ellps=GRS80 +no_defs"au_ll <- spTransform(pol, CRSobj = projection(lct))
# plot(pol_rs, border = "white", col = "gray", lwd = 2, axes = T)
#plot(clct)Correspondência entre células raster, geogrid e nomes dos binários
# referemcia da ordem dos tiles de cada binário
mat <- matrix(1:ncell(ref), ncol = ncol(ref))
guia_geogrid <- mat[nrow(mat):1,]
guia_geogrid## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 18 36 54 72 90 108 126 144 162 180 198 216 234
## [2,] 17 35 53 71 89 107 125 143 161 179 197 215 233
## [3,] 16 34 52 70 88 106 124 142 160 178 196 214 232
## [4,] 15 33 51 69 87 105 123 141 159 177 195 213 231
## [5,] 14 32 50 68 86 104 122 140 158 176 194 212 230
## [6,] 13 31 49 67 85 103 121 139 157 175 193 211 229
## [7,] 12 30 48 66 84 102 120 138 156 174 192 210 228
## [8,] 11 29 47 65 83 101 119 137 155 173 191 209 227
## [9,] 10 28 46 64 82 100 118 136 154 172 190 208 226
## [10,] 9 27 45 63 81 99 117 135 153 171 189 207 225
## [11,] 8 26 44 62 80 98 116 134 152 170 188 206 224
## [12,] 7 25 43 61 79 97 115 133 151 169 187 205 223
## [13,] 6 24 42 60 78 96 114 132 150 168 186 204 222
## [14,] 5 23 41 59 77 95 113 131 149 167 185 203 221
## [15,] 4 22 40 58 76 94 112 130 148 166 184 202 220
## [16,] 3 21 39 57 75 93 111 129 147 165 183 201 219
## [17,] 2 20 38 56 74 92 110 128 146 164 182 200 218
## [18,] 1 19 37 55 73 91 109 127 145 163 181 199 217
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 252 270 288 306 324 342 360 378 396 414 432
## [2,] 251 269 287 305 323 341 359 377 395 413 431
## [3,] 250 268 286 304 322 340 358 376 394 412 430
## [4,] 249 267 285 303 321 339 357 375 393 411 429
## [5,] 248 266 284 302 320 338 356 374 392 410 428
## [6,] 247 265 283 301 319 337 355 373 391 409 427
## [7,] 246 264 282 300 318 336 354 372 390 408 426
## [8,] 245 263 281 299 317 335 353 371 389 407 425
## [9,] 244 262 280 298 316 334 352 370 388 406 424
## [10,] 243 261 279 297 315 333 351 369 387 405 423
## [11,] 242 260 278 296 314 332 350 368 386 404 422
## [12,] 241 259 277 295 313 331 349 367 385 403 421
## [13,] 240 258 276 294 312 330 348 366 384 402 420
## [14,] 239 257 275 293 311 329 347 365 383 401 419
## [15,] 238 256 274 292 310 328 346 364 382 400 418
## [16,] 237 255 273 291 309 327 345 363 381 399 417
## [17,] 236 254 272 290 308 326 344 362 380 398 416
## [18,] 235 253 271 289 307 325 343 361 379 397 415
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 450 468 486 504 522 540 558 576 594 612 630
## [2,] 449 467 485 503 521 539 557 575 593 611 629
## [3,] 448 466 484 502 520 538 556 574 592 610 628
## [4,] 447 465 483 501 519 537 555 573 591 609 627
## [5,] 446 464 482 500 518 536 554 572 590 608 626
## [6,] 445 463 481 499 517 535 553 571 589 607 625
## [7,] 444 462 480 498 516 534 552 570 588 606 624
## [8,] 443 461 479 497 515 533 551 569 587 605 623
## [9,] 442 460 478 496 514 532 550 568 586 604 622
## [10,] 441 459 477 495 513 531 549 567 585 603 621
## [11,] 440 458 476 494 512 530 548 566 584 602 620
## [12,] 439 457 475 493 511 529 547 565 583 601 619
## [13,] 438 456 474 492 510 528 546 564 582 600 618
## [14,] 437 455 473 491 509 527 545 563 581 599 617
## [15,] 436 454 472 490 508 526 544 562 580 598 616
## [16,] 435 453 471 489 507 525 543 561 579 597 615
## [17,] 434 452 470 488 506 524 542 560 578 596 614
## [18,] 433 451 469 487 505 523 541 559 577 595 613
## [,36]
## [1,] 648
## [2,] 647
## [3,] 646
## [4,] 645
## [5,] 644
## [6,] 643
## [7,] 642
## [8,] 641
## [9,] 640
## [10,] 639
## [11,] 638
## [12,] 637
## [13,] 636
## [14,] 635
## [15,] 634
## [16,] 633
## [17,] 632
## [18,] 631# raster do geogrid
r_geogrid <- raster(ref)
r_geogrid[] <- guia_geogrid
plot(r_geogrid)
grid(nx=36,ny=18)
text(r_geogrid, cex = 0.5)
cels_correspondence <- data.frame(cel_raster = 1:ncell(r_geogrid),
cell_geog = r_geogrid[])
# componente direita do nome dos binarios
y_ini <- seq(from = 1, by = 1200, length.out = 18)
y_fim <- y_ini + 1200 -1
y_ini_char <- formatC(y_ini, digits = 0, width = 5, flag = '0', format = 'f')
y_fim_char <- formatC(y_fim, digits = 0, width = 5, flag = '0', format = 'f')
y <- paste0(y_ini_char, "-",y_fim_char)
# componente esquerda do nome dos binarios
x_ini <- seq(from = 1, by = 1200, length.out = 36)
x_fim <- x_ini + 1200 -1
x_ini_char <- formatC(x_ini, digits = 0, width = 5, flag = '0', format = 'f')
x_fim_char <- formatC(x_fim, digits = 0, width = 5, flag = '0', format = 'f')
x <- paste0(x_ini_char, "-",x_fim_char)
nomes_binarios <- expand.grid(x, y)
nomes_binarios <- do.call("paste", c(nomes_binarios, sep = "."))
nomes_binarios <- paste0(x, ".", y)
head(cbind(nomes_binarios), 20)## nomes_binarios
## [1,] "00001-01200.00001-01200"
## [2,] "01201-02400.01201-02400"
## [3,] "02401-03600.02401-03600"
## [4,] "03601-04800.03601-04800"
## [5,] "04801-06000.04801-06000"
## [6,] "06001-07200.06001-07200"
## [7,] "07201-08400.07201-08400"
## [8,] "08401-09600.08401-09600"
## [9,] "09601-10800.09601-10800"
## [10,] "10801-12000.10801-12000"
## [11,] "12001-13200.12001-13200"
## [12,] "13201-14400.13201-14400"
## [13,] "14401-15600.14401-15600"
## [14,] "15601-16800.15601-16800"
## [15,] "16801-18000.16801-18000"
## [16,] "18001-19200.18001-19200"
## [17,] "19201-20400.19201-20400"
## [18,] "20401-21600.20401-21600"
## [19,] "21601-22800.00001-01200"
## [20,] "22801-24000.01201-02400"tail(cbind(nomes_binarios), 20)## nomes_binarios
## [17,] "19201-20400.19201-20400"
## [18,] "20401-21600.20401-21600"
## [19,] "21601-22800.00001-01200"
## [20,] "22801-24000.01201-02400"
## [21,] "24001-25200.02401-03600"
## [22,] "25201-26400.03601-04800"
## [23,] "26401-27600.04801-06000"
## [24,] "27601-28800.06001-07200"
## [25,] "28801-30000.07201-08400"
## [26,] "30001-31200.08401-09600"
## [27,] "31201-32400.09601-10800"
## [28,] "32401-33600.10801-12000"
## [29,] "33601-34800.12001-13200"
## [30,] "34801-36000.13201-14400"
## [31,] "36001-37200.14401-15600"
## [32,] "37201-38400.15601-16800"
## [33,] "38401-39600.16801-18000"
## [34,] "39601-40800.18001-19200"
## [35,] "40801-42000.19201-20400"
## [36,] "42001-43200.20401-21600"files_cells_corresp <- data.frame(nomes_binarios, cels_correspondence)
select_files <- dplyr::filter(files_cells_corresp, cel_raster %in% cels)
select_files| nomes_binarios | cel_raster | cell_geog |
|---|---|---|
| 14401-15600.14401-15600 | 409 | 223 |
| 15601-16800.15601-16800 | 410 | 241 |
| 14401-15600.14401-15600 | 445 | 222 |
Dados binários
- descobrir tile e nome do arquivo que contem o dominio espacial do RS.
- com o nome do arquivo importar no R/NCL para visualizar as categorias (20 categorias), mas no RS nao deve conter as categorias 18-20 (de tundra) que foram adicionadas em relação ao mapa de LUC do MODIS (17 categorias).
- descobrir o formato do binário (quantos bytes (1,…,4)?, inteiro, big-endian byte order, 1200*1200)
x <- readBin(con = "../../atividade1/data/input_data_wrf/modis_landuse_20class_30s/33601-34800.09601-10800",
what = "integer",
size = 1,
endian = "big",
n = 1200*1200)
x[1:100]## [1] 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
## [24] 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
## [47] 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
## [70] 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
## [93] 17 17 17 17 17 17 17 17range(x)## [1] 1 17length(x)## [1] 1440000hist(x)
gerar um raster com áreas úmidas para o tile que contém o RS. Precisa saber a “extent” (domínio espacial do tile)
sobrescrever as células dos dados da FEPAM nas células correspondentes do binário do tile