Diff to HTML by rtfpessoa

Files changed (1) show
{orig → migrated_terra}/die_mod.R RENAMED
@@ -4,27 +4,40 @@
4
4
  ###################################################
5
5
  ### code chunk number 9: die.Rnw:95-96
6
6
  ###################################################
7
- library(rgdal)
8
7
 
9
8
 
10
9
  ###################################################
11
10
  ### code chunk number 14: die.Rnw:264-266
12
11
  ###################################################
13
12
  # From PROJ 6.0.0, EPSG data stored in an SQLite database without proj4 strings
14
- EPSG <- try(make_EPSG())
15
- if (class(EPSG) != "try-error") EPSG[grep("^# ED50$", EPSG$note),]
13
+
14
+ pths <- sf::sf_proj_search_paths() # terra without access to search path
15
+ library(RSQLite)
16
+ db <- dbConnect(SQLite(), dbname=file.path(pths[length(pths)], "proj.db"))
17
+ GCRSi <- dbGetQuery(db, "SELECT auth_name, code, name FROM geodetic_crs WHERE typeof(code) = 'integer'")
18
+ GCRSi$code <- as.character(GCRSi$code)
19
+ GCRSt <- dbGetQuery(db, "SELECT auth_name, code, name FROM geodetic_crs WHERE typeof(code) = 'text'")
20
+ PCRSi <- dbGetQuery(db, "SELECT auth_name, code, name FROM projected_crs WHERE typeof(code) = 'integer'")
21
+ PCRSi$code <- as.character(PCRSi$code)
22
+ PCRSt <- dbGetQuery(db, "SELECT auth_name, code, name FROM projected_crs WHERE typeof(code) = 'text'")
23
+ dbDisconnect(db)
24
+ GPCRS <- rbind(GCRSi, GCRSt, PCRSi, PCRSt)
25
+ GPCRS[grep("^ED50$", GPCRS$name),]
16
26
 
17
27
 
18
28
  ###################################################
19
29
  ### code chunk number 15: die.Rnw:306-307
20
30
  ###################################################
21
- CRS("+init=epsg:4230")
31
+ library(terra)
32
+ cat(crs("EPSG:4230"), "\n")
22
33
 
23
34
 
24
35
  ###################################################
25
36
  ### code chunk number 16: die.Rnw:330-332
26
37
  ###################################################
27
- ED50 <- CRS("+init=epsg:4230 +towgs84=-87,-96,-120,0,0,0,0")
38
+ library(sp)
39
+ set_evolution_status(2L)
40
+ ED50 <- CRS("EPSG:4230")
28
41
  ED50
29
42
 
30
43
 
@@ -34,43 +47,40 @@
34
47
  IJ.east <- as(char2dms("4d31\'00\"E"), "numeric")
35
48
  IJ.north <- as(char2dms("52d28\'00\"N"), "numeric")
36
49
  IJ.ED50 <- SpatialPoints(cbind(x=IJ.east, y=IJ.north), proj4string=ED50)
37
- res <- spTransform(IJ.ED50, CRS("+proj=longlat +datum=WGS84"))
50
+ res <- as(sf::st_transform(sf::st_as_sf(IJ.ED50), sf::st_crs(CRS("EPSG:4326"))), "Spatial")
38
51
  x <- as(dd2dms(coordinates(res)[1]), "character")
39
52
  y <- as(dd2dms(coordinates(res)[2], TRUE), "character")
40
53
  cat(x, y, "\n")
41
54
  spDistsN1(coordinates(IJ.ED50), coordinates(res), longlat=TRUE)*1000
42
- library(maptools)
43
- gzAzimuth(coordinates(IJ.ED50), coordinates(res))
55
+ p = sf::st_sfc(sf::st_point(coordinates(IJ.ED50)), sf::st_point(coordinates(res)), crs = 4326)
56
+ lwgeom::st_geod_azimuth(p) |> units::set_units("degree")
44
57
 
45
58
 
46
59
  ###################################################
47
60
  ### code chunk number 19: die.Rnw:430-434
48
61
  ###################################################
49
- proj4string(IJ.ED50) <- CRS("+init=epsg:4230")
50
- res <- spTransform(IJ.ED50, CRS("+proj=longlat +datum=WGS84"))
62
+ proj4string(IJ.ED50) <- CRS("EPSG:4230")
63
+ res <- as(sf::st_transform(sf::st_as_sf(IJ.ED50), sf::st_crs(CRS("EPSG:4326"))), "Spatial")
51
64
  spDistsN1(coordinates(IJ.ED50), coordinates(res), longlat=TRUE)*1000
52
- gzAzimuth(coordinates(IJ.ED50), coordinates(res))
65
+ p = sf::st_sfc(sf::st_point(coordinates(IJ.ED50)), sf::st_point(coordinates(res)), crs = 4326)
66
+ lwgeom::st_geod_azimuth(p) |> units::set_units("degree")
53
67
 
54
68
 
55
69
  ###################################################
56
70
  ### code chunk number 21: die.Rnw:449-450
57
71
  ###################################################
58
- if (class(EPSG) != "try-error") EPSG[grep("Atlas", EPSG$note), 1:2]
72
+ GPCRS[grep("US National Atlas", GPCRS$name),]
59
73
 
60
74
 
61
75
  ###################################################
62
76
  ### code chunk number 22: die.Rnw:452-453
63
77
  ###################################################
64
- CRS("+init=epsg:2163")
78
+ cat(crs("EPSG:2163"), "\n")
65
79
 
66
80
 
67
81
  ###################################################
68
82
  ### code chunk number 24: die.Rnw:475-479
69
83
  ###################################################
70
- proj <- projInfo("proj")
71
- proj[proj$name == "laea",]
72
- ellps <- projInfo("ellps")
73
- ellps[grep("a=6370997", ellps$major),]
74
84
 
75
85
 
76
86
  ###################################################
@@ -99,7 +109,7 @@
99
109
  ###################################################
100
110
  ### code chunk number 28: die.Rnw:672-673
101
111
  ###################################################
102
- head(ogrDrivers(), n=10)
112
+ head(gdal(drivers=TRUE) |> subset(subset=type=="vector"), n=10)
103
113
 
104
114
 
105
115
  ###################################################
@@ -112,15 +122,15 @@
112
122
  ###################################################
113
123
  ### code chunk number 34: die.Rnw:760-761
114
124
  ###################################################
115
- ogrInfo(".", "scot")
116
125
 
117
126
 
118
127
  ###################################################
119
128
  ### code chunk number 37: die.Rnw:791-794
120
129
  ###################################################
121
- scot_LL <- readOGR(dsn="scot.shp", layer="scot", integer64="allow.loss")
130
+ library(raster)
131
+ scot_LL <- as(vect("scot.shp"), "Spatial")
122
132
  proj4string(scot_LL)
123
- proj4string(scot_LL) <- CRS("+proj=longlat +ellps=WGS84")
133
+ proj4string(scot_LL) <- CRS("EPSG:4326")
124
134
 
125
135
 
126
136
  ###################################################
@@ -137,9 +147,7 @@
137
147
  ID_D <- match(scot_LL$ID, scot_dat$District)
138
148
  scot_dat1 <- scot_dat[ID_D,]
139
149
  row.names(scot_dat1) <- row.names(scot_LL)
140
- library(maptools)
141
- scot_LLa <- spCbind(scot_LL, scot_dat1)
142
- all.equal(scot_LLa$ID, scot_LLa$District)
150
+ scot_LLa <- as(merge(vect(scot_LL), scot_dat1, by.x="ID", by.y="District"), "Spatial")
143
151
  names(scot_LLa)
144
152
 
145
153
 
@@ -150,14 +158,13 @@
150
158
  O <- scot_LLa$Observed
151
159
  E <- scot_LLa$Expected
152
160
  scot_LLa$SMR <- probmap(O, E)$relRisk/100
153
- library(DCluster)
154
- scot_LLa$smth <- empbaysmooth(O, E)$smthrr
161
+ scot_LLa$smth <- DCluster::empbaysmooth(O, E)$smthrr
155
162
 
156
163
 
157
164
  ###################################################
158
165
  ### code chunk number 42: die.Rnw:863-864
159
166
  ###################################################
160
- scot_BNG <- spTransform(scot_LLa, CRS("+init=epsg:27700"))
167
+ scot_BNG <- as(project(vect(scot_LLa), "EPSG:27700"), "Spatial")
161
168
 
162
169
 
163
170
  ###################################################
@@ -177,8 +184,7 @@
177
184
  ###################################################
178
185
  ### code chunk number 45: die.Rnw:907-909
179
186
  ###################################################
180
- drv <- "ESRI Shapefile"
181
- writeOGR(scot_BNG, dsn=".", layer="scot_BNG", driver=drv, overwrite_layer=TRUE)
187
+ writeVector(vect(scot_BNG), "scot_BNG.shp", overwrite=TRUE)
182
188
 
183
189
 
184
190
  ###################################################
@@ -221,7 +227,7 @@
221
227
  ###################################################
222
228
  ### code chunk number 53: die.Rnw:961-962
223
229
  ###################################################
224
- Fires <- readOGR("fires_120104.shp", "fires_120104")
230
+ Fires <- as(vect("fires_120104.shp"), "Spatial")
225
231
 
226
232
 
227
233
  ###################################################
@@ -236,12 +242,11 @@
236
242
  x <- c(-15, -15, 38, 38, -15)
237
243
  y <- c(28, 62, 62, 28, 28)
238
244
  crds <- cbind(x=x, y=y)
239
- bb <- SpatialPolygons(list(Polygons(list(Polygon(coords=crds)), "1")))
240
- library(maptools)
241
- data(wrld_simpl)
242
- proj4string(bb) <- CRS(proj4string(wrld_simpl))
243
- library(rgeos)
244
- slbb <- gIntersection(bb, as(wrld_simpl, "SpatialLines"))
245
+ bb <- vect(SpatialPolygons(list(Polygons(list(Polygon(coords=crds)), "1"))))
246
+ data(world, package="spData")
247
+ world <- as(world, "Spatial")
248
+ crs(bb) <- wkt(world)
249
+ slbb <- as(crop(vect(as(world, "SpatialLines")), bb), "Spatial")
245
250
  spl <- list("sp.lines", slbb, lwd=0.7, col="khaki4")
246
251
 
247
252
 
@@ -281,13 +286,13 @@
281
286
  ###################################################
282
287
  ### code chunk number 63: die.Rnw:1107-1108
283
288
  ###################################################
284
- getinfo.shape("scot_BNG.shp")
285
289
 
286
290
 
287
291
  ###################################################
288
292
  ### code chunk number 66: die.Rnw:1178-1181
289
293
  ###################################################
290
- auck_el1 <- readGDAL("70042108.tif")
294
+ auck_el1 <- as(raster(rast("70042108.tif")), "SpatialGridDataFrame")
295
+ names(auck_el1) <- "band1"
291
296
  summary(auck_el1)
292
297
  is.na(auck_el1$band1) <- auck_el1$band1 <= 0 | auck_el1$band1 > 1e+4
293
298
 
@@ -295,13 +300,6 @@
295
300
  ###################################################
296
301
  ### code chunk number 69: die.Rnw:1204-1211
297
302
  ###################################################
298
- x <- GDAL.open("70042108.tif")
299
- xx <- getDriver(x)
300
- #xx #do not show pointer
301
- getDriverLongName(xx)
302
- #x #do not show pointer
303
- dim(x)
304
- GDAL.close(x)
305
303
 
306
304
 
307
305
  ###################################################
@@ -312,10 +310,10 @@
312
310
  pal
313
311
  length(pal) == length(brks)-1
314
312
  auck_el1$band1 <- findInterval(auck_el1$band1, vec=brks, all.inside=TRUE)-1
315
- writeGDAL(auck_el1, "demIndex.tif", drivername="GTiff", type="Byte", colorTable=list(pal), mvFlag=length(brks)-1)
316
- Gi <- GDALinfo("demIndex.tif", returnColorTable=TRUE)
317
- CT <- attr(Gi, "ColorTable")[[1]]
318
- CT[CT > "#000000"]
313
+ o <- rast(auck_el1)
314
+ coltab(o, layer=1) <- data.frame(value=0:10, col=pal)
315
+ writeRaster(o, "demIndex.tif", filetype="GTiff", datatype="INT1U", NAflag=length(brks)-1, overwrite=TRUE)
316
+ head(coltab(rast("demIndex.tif"))[[1]], n=10)
319
317
 
320
318
 
321
319
  ###################################################
@@ -324,7 +322,7 @@
324
322
  data(meuse.grid)
325
323
  coordinates(meuse.grid) <- c("x", "y")
326
324
  gridded(meuse.grid) <- TRUE
327
- slot(meuse.grid, "proj4string") <- CRS("+init=epsg:28992")
325
+ slot(meuse.grid, "proj4string") <- CRS("EPSG:28992")
328
326
  data(meuse)
329
327
  coordinates(meuse) <- c("x", "y")
330
328
  slot(meuse, "proj4string") <- slot(meuse.grid, "proj4string")
@@ -341,13 +339,13 @@
341
339
  ### code chunk number 78: die.Rnw:1312-1314
342
340
  ###################################################
343
341
  summary(log_zinc)
344
- writeGDAL(log_zinc, fname="log_zinc.tif", drivername="GTiff", type="Float32", options="INTERLEAVE=PIXEL")
342
+ writeRaster(rast(log_zinc), filename="log_zinc.tif", filetype="GTiff", datatype="FLT4S", gdal="INTERLEAVE=PIXEL", overwrite=TRUE)
345
343
 
346
344
 
347
345
  ###################################################
348
346
  ### code chunk number 79: die.Rnw:1316-1317
349
347
  ###################################################
350
- GDALinfo("log_zinc.tif")
348
+ rast("log_zinc.tif")
351
349
 
352
350
 
353
351
  ###################################################
@@ -355,18 +353,16 @@
355
353
  ###################################################
356
354
  Soil <- meuse.grid["soil"]
357
355
  table(Soil$soil)
358
- Soil$soil <- as.integer(Soil$soil)-1
359
356
  Cn <- c("Rd10A", "Rd90C/VII", "Bkd26/VII")
360
- writeGDAL(Soil, "Soil.tif", drivername="GTiff", type="Byte", catNames=list(Cn), mvFlag=length(Cn))
361
- Gi <- GDALinfo("Soil.tif", returnCategoryNames=TRUE)
362
- attr(Gi, "CATlist")[[1]]
363
- summary(readGDAL("Soil.tif"))
357
+ Soil$soil <- factor(Soil$soil, levels=1:3, labels=Cn)
358
+ writeRaster(rast(Soil), filename="Soil.tif", filetype="GTiff", datatype="INT1U", overwrite=TRUE)
359
+ summary(rast("Soil.tif"))
364
360
 
365
361
 
366
362
  ###################################################
367
363
  ### code chunk number 84: die.Rnw:1368-1369
368
364
  ###################################################
369
- head(gdalDrivers(), n=10)
365
+ head(gdal(drivers=TRUE) |> subset(subset=type=="raster"), n=10)
370
366
 
371
367
 
372
368
  ###################################################
@@ -387,19 +383,19 @@
387
383
  ###################################################
388
384
  ### code chunk number 94: die.Rnw:1425-1426
389
385
  ###################################################
390
- osm <- readGDAL("osm_bergen_120105.tif")
386
+ osm123 <- rast("osm_bergen_120105.tif")
391
387
 
392
388
 
393
389
  ###################################################
394
390
  ### code chunk number 96: die.Rnw:1431-1432
395
391
  ###################################################
396
- summary(osm)
392
+ summary(osm123)
397
393
 
398
394
 
399
395
  ###################################################
400
396
  ### code chunk number 97: die.Rnw:1441-1450
401
397
  ###################################################
402
- image(osm, red=1, green=2, blue=3)
398
+ plot(osm123)
403
399
 
404
400
 
405
401
  ###################################################
@@ -424,7 +420,7 @@
424
420
  cs <- dBB/DIM12
425
421
  cc <- c(BB[1,2] + cs[1]/2, BB[1,1] + cs[2]/2)
426
422
  GT <- GridTopology(cc, cs, DIM12)
427
- p4s <- CRS("+proj=longlat +datum=WGS84")
423
+ p4s <- CRS("EPSG:4326")
428
424
  SG_myMap <- SpatialGridDataFrame(GT, proj4string=p4s, data=data.frame(r=c(t(myMap$myTile[,,1]))*255, g=c(t(myMap$myTile[,,2]))*255, b=c(t(myMap$myTile[,,3]))*255))
429
425
 
430
426
 
@@ -437,8 +433,8 @@
437
433
  ###################################################
438
434
  ### code chunk number 102: die.Rnw:1535-1537
439
435
  ###################################################
440
- library(osmar)
441
- load("osmar.RData")
436
+ #library(osmar)
437
+ #load("osmar.RData")
442
438
 
443
439
 
444
440
  ###################################################
@@ -453,76 +449,78 @@
453
449
  ###################################################
454
450
  ### code chunk number 104: die.Rnw:1546-1547
455
451
  ###################################################
456
- torget1 <- as_sp(torget, "lines")
452
+ #torget1 <- as_sp(torget, "lines")
457
453
 
458
454
 
459
455
  ###################################################
460
456
  ### code chunk number 105: die.Rnw:1548-1549
461
457
  ###################################################
462
- sort(table(torget1$user), decreasing=TRUE)[1:3]
458
+ #sort(table(torget1$user), decreasing=TRUE)[1:3]
463
459
 
464
460
 
465
461
  ###################################################
466
462
  ### code chunk number 107: die.Rnw:1572-1576
467
463
  ###################################################
468
- bybane <- find(torget, way(tags(k == "light_rail")))
469
- bybane <- find_down(torget, way(bybane))
470
- bybane <- subset(torget, ids=bybane)
471
- bybane <- as_sp(bybane, "lines")
464
+ #bybane <- find(torget, way(tags(k == "light_rail")))
465
+ #bybane <- find_down(torget, way(bybane))
466
+ #bybane <- subset(torget, ids=bybane)
467
+ #bybane <- as_sp(bybane, "lines")
472
468
 
473
469
 
474
470
  ###################################################
475
471
  ### code chunk number 108: die.Rnw:1583-1596
476
472
  ###################################################
477
- image(SG_myMap, red=1, green=2, blue=3)
478
- plot(torget1, add=TRUE)
479
- plot(bybane, add=TRUE, lwd=5, col="orange2")
480
- plot(0:1, 0:1, type = "n", axes = FALSE, asp=1)
481
- rasterImage(myMap1[[4]], 0, 0, 1, 1)
473
+ #image(SG_myMap, red=1, green=2, blue=3)
474
+ #plot(torget1, add=TRUE)
475
+ #plot(bybane, add=TRUE, lwd=5, col="orange2")
476
+ #plot(0:1, 0:1, type = "n", axes = FALSE, asp=1)
477
+ #rasterImage(myMap1[[4]], 0, 0, 1, 1)
482
478
 
483
479
 
484
480
  ###################################################
485
481
  ### code chunk number 110: die.Rnw:1614-1615
486
482
  ###################################################
487
- writeOGR(Fires[,c("gml_id", "FireDate", "Area_HA")], dsn="fires.kml", layer="fires", driver="KML", overwrite_layer=TRUE)
483
+ writeVector(vect(Fires[,c("gml_id", "FireDate", "Area_HA")]), filename="fires.kml", layer="fires", filetype="KML", options=NULL, overwrite=TRUE)
488
484
 
489
485
 
490
486
  ###################################################
491
487
  ### code chunk number 111: die.Rnw:1673-1679
492
488
  ###################################################
493
- library(maptools)
494
489
  grd <- as(meuse.grid, "SpatialPolygons")
495
490
  proj4string(grd) <- CRS(proj4string(meuse))
496
- grd.union <- unionSpatialPolygons(grd, rep("x", length(slot(grd, "polygons"))))
497
- ll <- CRS("+proj=longlat +datum=WGS84")
498
- grd.union.ll <- spTransform(grd.union, ll)
491
+ grd.union <- aggregate(vect(grd))
492
+ grd.union.ll <- project(grd.union, "EPSG:4326")
499
493
 
500
494
 
501
495
  ###################################################
502
496
  ### code chunk number 112: die.Rnw:1703-1707
503
497
  ###################################################
504
- llGRD <- GE_SpatialGrid(grd.union.ll)
505
- llGRD_in <- over(llGRD$SG, grd.union.ll)
506
- llSGDF <- SpatialGridDataFrame(grid=slot(llGRD$SG, "grid"), proj4string=CRS(proj4string(llGRD$SG)), data=data.frame(in0=llGRD_in))
498
+ llGRD <- rast(ext(grd.union.ll), ncol=600, nrow=600)
499
+ values(llGRD) <- 1
500
+ llGRD_in <- mask(llGRD, grd.union.ll)
501
+ llSGDF <- as(raster(llGRD_in), "SpatialGridDataFrame")
507
502
  llSPix <- as(llSGDF, "SpatialPixelsDataFrame")
508
503
 
509
504
 
510
505
  ###################################################
511
506
  ### code chunk number 113: die.Rnw:1722-1724
512
507
  ###################################################
513
- meuse_ll <- spTransform(meuse, CRS("+proj=longlat +datum=WGS84"))
508
+ meuse_ll <- as(project(vect(meuse), CRS("EPSG:4326")), "Spatial")
514
509
  llSPix$pred <- gstat::idw(log(zinc)~1, meuse_ll, llSPix)$var1.pred
515
510
 
516
511
 
517
512
  ###################################################
518
513
  ### code chunk number 114: die.Rnw:1744-1749
519
514
  ###################################################
520
- png(file="zinc_IDW.png", width=llGRD$width, height=llGRD$height, bg="transparent")
521
- par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
522
515
  image(llSPix, "pred", col=bpy.colors(20))
523
- dev.off()
524
- kmlOverlay(llGRD, "zinc_IDW.kml", "zinc_IDW.png")
516
+
525
517
 
526
518
  ###################################################
527
- source("die_snow.R", echo=TRUE)
519
+ #source("die_snow.R", echo=TRUE)
520
+
521
+ (sI <- sessionInfo()) # check: no sp?
522
+
523
+ "rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
524
+ "rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
525
+ "maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
528
526