Diff to HTML by rtfpessoa

Files changed (1) show
{orig → migrated}/die_mod.R RENAMED
@@ -4,27 +4,37 @@
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),]
16
13
 
14
+ library(sf)
15
+ pths <- sf_proj_search_paths()
16
+ library(RSQLite)
17
+ db <- dbConnect(SQLite(), dbname=file.path(pths[length(pths)], "proj.db"))
18
+ GCRSi <- dbGetQuery(db, "SELECT auth_name, code, name FROM geodetic_crs WHERE typeof(code) = 'integer'")
19
+ GCRSi$code <- as.character(GCRSi$code)
20
+ GCRSt <- dbGetQuery(db, "SELECT auth_name, code, name FROM geodetic_crs WHERE typeof(code) = 'text'")
21
+ PCRSi <- dbGetQuery(db, "SELECT auth_name, code, name FROM projected_crs WHERE typeof(code) = 'integer'")
22
+ PCRSi$code <- as.character(PCRSi$code)
23
+ PCRSt <- dbGetQuery(db, "SELECT auth_name, code, name FROM projected_crs WHERE typeof(code) = 'text'")
24
+ dbDisconnect(db)
25
+ GPCRS <- rbind(GCRSi, GCRSt, PCRSi, PCRSt)
26
+ GPCRS[grep("^ED50$", GPCRS$name),]
17
27
 
18
28
  ###################################################
19
29
  ### code chunk number 15: die.Rnw:306-307
20
30
  ###################################################
21
- CRS("+init=epsg:4230")
22
-
31
+ library(sp)
32
+ CRS("EPSG:4230")
23
33
 
24
34
  ###################################################
25
35
  ### code chunk number 16: die.Rnw:330-332
26
36
  ###################################################
27
- ED50 <- CRS("+init=epsg:4230 +towgs84=-87,-96,-120,0,0,0,0")
37
+ ED50 <- CRS("EPSG:4230")
28
38
  ED50
29
39
 
30
40
 
@@ -34,43 +44,38 @@
34
44
  IJ.east <- as(char2dms("4d31\'00\"E"), "numeric")
35
45
  IJ.north <- as(char2dms("52d28\'00\"N"), "numeric")
36
46
  IJ.ED50 <- SpatialPoints(cbind(x=IJ.east, y=IJ.north), proj4string=ED50)
37
- res <- spTransform(IJ.ED50, CRS("+proj=longlat +datum=WGS84"))
47
+ res <- as(st_transform(st_as_sf(IJ.ED50), st_crs(CRS("EPSG:4326"))), "Spatial")
38
48
  x <- as(dd2dms(coordinates(res)[1]), "character")
39
49
  y <- as(dd2dms(coordinates(res)[2], TRUE), "character")
40
50
  cat(x, y, "\n")
41
51
  spDistsN1(coordinates(IJ.ED50), coordinates(res), longlat=TRUE)*1000
42
- library(maptools)
43
- gzAzimuth(coordinates(IJ.ED50), coordinates(res))
44
-
52
+ p = st_sfc(st_point(coordinates(IJ.ED50)), st_point(coordinates(res)), crs = 4326)
53
+ library(lwgeom)
54
+ st_geod_azimuth(p) |> units::set_units("degree")
45
55
 
46
56
  ###################################################
47
57
  ### code chunk number 19: die.Rnw:430-434
48
58
  ###################################################
49
- proj4string(IJ.ED50) <- CRS("+init=epsg:4230")
50
- res <- spTransform(IJ.ED50, CRS("+proj=longlat +datum=WGS84"))
59
+ proj4string(IJ.ED50) <- CRS("EPSG:4230")
60
+ res <- as(st_transform(st_as_sf(IJ.ED50), st_crs(CRS("EPSG:4326"))), "Spatial")
51
61
  spDistsN1(coordinates(IJ.ED50), coordinates(res), longlat=TRUE)*1000
52
- gzAzimuth(coordinates(IJ.ED50), coordinates(res))
53
-
62
+ p = st_sfc(st_point(coordinates(IJ.ED50)), st_point(coordinates(res)), crs = 4326)
63
+ st_geod_azimuth(p) |> units::set_units("degree")
54
64
 
55
65
  ###################################################
56
66
  ### code chunk number 21: die.Rnw:449-450
57
67
  ###################################################
58
- if (class(EPSG) != "try-error") EPSG[grep("Atlas", EPSG$note), 1:2]
59
-
68
+ GPCRS[grep("US National Atlas", GPCRS$name),]
60
69
 
61
70
  ###################################################
62
71
  ### code chunk number 22: die.Rnw:452-453
63
72
  ###################################################
64
- CRS("+init=epsg:2163")
73
+ CRS("EPSG:2163")
65
74
 
66
75
 
67
76
  ###################################################
68
77
  ### code chunk number 24: die.Rnw:475-479
69
78
  ###################################################
70
- proj <- projInfo("proj")
71
- proj[proj$name == "laea",]
72
- ellps <- projInfo("ellps")
73
- ellps[grep("a=6370997", ellps$major),]
74
79
 
75
80
 
76
81
  ###################################################
@@ -99,7 +104,7 @@
99
104
  ###################################################
100
105
  ### code chunk number 28: die.Rnw:672-673
101
106
  ###################################################
102
- head(ogrDrivers(), n=10)
107
+ head(st_drivers(), n=10)
103
108
 
104
109
 
105
110
  ###################################################
@@ -112,15 +117,14 @@
112
117
  ###################################################
113
118
  ### code chunk number 34: die.Rnw:760-761
114
119
  ###################################################
115
- ogrInfo(".", "scot")
116
120
 
117
121
 
118
122
  ###################################################
119
123
  ### code chunk number 37: die.Rnw:791-794
120
124
  ###################################################
121
- scot_LL <- readOGR(dsn="scot.shp", layer="scot", integer64="allow.loss")
125
+ scot_LL <- as(st_read("scot.shp"), "Spatial")
122
126
  proj4string(scot_LL)
123
- proj4string(scot_LL) <- CRS("+proj=longlat +ellps=WGS84")
127
+ proj4string(scot_LL) <- CRS("EPSG:4326")
124
128
 
125
129
 
126
130
  ###################################################
@@ -137,16 +141,14 @@
137
141
  ID_D <- match(scot_LL$ID, scot_dat$District)
138
142
  scot_dat1 <- scot_dat[ID_D,]
139
143
  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)
144
+ scot_LLa <- as(merge(st_as_sf(scot_LL), scot_dat1, by.x="ID", by.y="District"), "Spatial")
143
145
  names(scot_LLa)
144
146
 
145
147
 
146
148
  ###################################################
147
149
  ### code chunk number 41: die.Rnw:846-852
148
150
  ###################################################
149
- library(spdep)
151
+ library(spdep) # drags in raster via spData
150
152
  O <- scot_LLa$Observed
151
153
  E <- scot_LLa$Expected
152
154
  scot_LLa$SMR <- probmap(O, E)$relRisk/100
@@ -157,7 +159,7 @@
157
159
  ###################################################
158
160
  ### code chunk number 42: die.Rnw:863-864
159
161
  ###################################################
160
- scot_BNG <- spTransform(scot_LLa, CRS("+init=epsg:27700"))
162
+ scot_BNG <- as(st_transform(st_as_sf(scot_LLa), st_crs(CRS("EPSG:27700"))), "Spatial")
161
163
 
162
164
 
163
165
  ###################################################
@@ -177,8 +179,7 @@
177
179
  ###################################################
178
180
  ### code chunk number 45: die.Rnw:907-909
179
181
  ###################################################
180
- drv <- "ESRI Shapefile"
181
- writeOGR(scot_BNG, dsn=".", layer="scot_BNG", driver=drv, overwrite_layer=TRUE)
182
+ st_write(st_as_sf(scot_BNG), "scot_BNG.shp", delete_layer=TRUE)
182
183
 
183
184
 
184
185
  ###################################################
@@ -221,8 +222,7 @@
221
222
  ###################################################
222
223
  ### code chunk number 53: die.Rnw:961-962
223
224
  ###################################################
224
- Fires <- readOGR("fires_120104.shp", "fires_120104")
225
-
225
+ Fires <- as(st_read("fires_120104.shp"), "Spatial")
226
226
 
227
227
  ###################################################
228
228
  ### code chunk number 55: die.Rnw:967-968
@@ -236,12 +236,13 @@
236
236
  x <- c(-15, -15, 38, 38, -15)
237
237
  y <- c(28, 62, 62, 28, 28)
238
238
  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"))
239
+ bb <- st_as_sf(SpatialPolygons(list(Polygons(list(Polygon(coords=crds)), "1"))))
240
+ data(world, package="spData") # spData drags in raster
241
+ st_crs(bb) <- st_crs(world)
242
+ s2_status <- sf_use_s2()
243
+ sf_use_s2(FALSE)
244
+ slbb <- as(st_intersection(bb, st_cast(world, "MULTILINESTRING")), "Spatial")
245
+ sf_use_s2(s2_status)
245
246
  spl <- list("sp.lines", slbb, lwd=0.7, col="khaki4")
246
247
 
247
248
 
@@ -268,7 +269,6 @@
268
269
  GR_Fires <- Fires1[Fires1$Country == "GR",]
269
270
  #writeOGR(GR_Fires, "EFFIS.gpx", "waypoints", driver="GPX", dataset_options="GPX_USE_EXTENSIONS=YES", overwrite_layer=TRUE, delete_dsn=TRUE)
270
271
  #st_write(st_as_sf(GR_Fires), dsn="EFFIS.gpx", layer="waypoints", driver="GPX", dataset_options="GPX_USE_EXTENSIONS=YES", delete_layer=TRUE, delete_dsn=TRUE)
271
- # GDAL >= 3.6 GPX driver changes
272
272
 
273
273
  ###################################################
274
274
  ### code chunk number 61: die.Rnw:1071-1073
@@ -277,17 +277,19 @@
277
277
  #GR <- as(st_read("EFFIS.gpx", "waypoints"), "Spatial")
278
278
  #GR[1,c(5,24:28)]
279
279
  #base::print(GR[1,c(5,24:28)])
280
+ # skip GPX after GDAL > 3.6 driver change
280
281
 
281
282
  ###################################################
282
283
  ### code chunk number 63: die.Rnw:1107-1108
283
284
  ###################################################
284
- getinfo.shape("scot_BNG.shp")
285
285
 
286
286
 
287
287
  ###################################################
288
288
  ### code chunk number 66: die.Rnw:1178-1181
289
289
  ###################################################
290
- auck_el1 <- readGDAL("70042108.tif")
290
+ library(stars)
291
+ auck_el1 <- as(read_stars("70042108.tif"), "Spatial")
292
+ names(auck_el1) <- "band1"
291
293
  summary(auck_el1)
292
294
  is.na(auck_el1$band1) <- auck_el1$band1 <= 0 | auck_el1$band1 > 1e+4
293
295
 
@@ -295,27 +297,11 @@
295
297
  ###################################################
296
298
  ### code chunk number 69: die.Rnw:1204-1211
297
299
  ###################################################
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
300
 
306
301
 
307
302
  ###################################################
308
303
  ### code chunk number 74: die.Rnw:1263-1272
309
304
  ###################################################
310
- brks <- c(0,10,20,50,100,150,200,300,400,500,600,700)
311
- pal <- terrain.colors(11)
312
- pal
313
- length(pal) == length(brks)-1
314
- 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"]
319
305
 
320
306
 
321
307
  ###################################################
@@ -324,7 +310,7 @@
324
310
  data(meuse.grid)
325
311
  coordinates(meuse.grid) <- c("x", "y")
326
312
  gridded(meuse.grid) <- TRUE
327
- slot(meuse.grid, "proj4string") <- CRS("+init=epsg:28992")
313
+ slot(meuse.grid, "proj4string") <- CRS("EPSG:28992")
328
314
  data(meuse)
329
315
  coordinates(meuse) <- c("x", "y")
330
316
  slot(meuse, "proj4string") <- slot(meuse.grid, "proj4string")
@@ -341,13 +327,11 @@
341
327
  ### code chunk number 78: die.Rnw:1312-1314
342
328
  ###################################################
343
329
  summary(log_zinc)
344
- writeGDAL(log_zinc, fname="log_zinc.tif", drivername="GTiff", type="Float32", options="INTERLEAVE=PIXEL")
345
-
330
+ write_stars(st_as_stars(log_zinc), dsn="log_zinc.tif", type="Float32", options="INTERLEAVE=PIXEL")
346
331
 
347
332
  ###################################################
348
333
  ### code chunk number 79: die.Rnw:1316-1317
349
334
  ###################################################
350
- GDALinfo("log_zinc.tif")
351
335
 
352
336
 
353
337
  ###################################################
@@ -355,18 +339,15 @@
355
339
  ###################################################
356
340
  Soil <- meuse.grid["soil"]
357
341
  table(Soil$soil)
358
- Soil$soil <- as.integer(Soil$soil)-1
359
342
  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"))
343
+ Soil$soil <- factor(Soil$soil, levels=1:3, labels=Cn)
344
+ write_stars(st_as_stars(Soil), "Soil.tif", type="Byte")
364
345
 
365
346
 
366
347
  ###################################################
367
348
  ### code chunk number 84: die.Rnw:1368-1369
368
349
  ###################################################
369
- head(gdalDrivers(), n=10)
350
+ head(st_drivers("raster"), n=10)
370
351
 
371
352
 
372
353
  ###################################################
@@ -387,8 +368,8 @@
387
368
  ###################################################
388
369
  ### code chunk number 94: die.Rnw:1425-1426
389
370
  ###################################################
390
- osm <- readGDAL("osm_bergen_120105.tif")
391
-
371
+ osm123 <- read_stars("osm_bergen_120105.tif")
372
+ osm <- cbind(as(osm123[,,,1, drop=TRUE], "Spatial"), as(osm123[,,,2, drop=TRUE], "Spatial"), as(osm123[,,,3, drop=TRUE], "Spatial"))
392
373
 
393
374
  ###################################################
394
375
  ### code chunk number 96: die.Rnw:1431-1432
@@ -424,7 +405,7 @@
424
405
  cs <- dBB/DIM12
425
406
  cc <- c(BB[1,2] + cs[1]/2, BB[1,1] + cs[2]/2)
426
407
  GT <- GridTopology(cc, cs, DIM12)
427
- p4s <- CRS("+proj=longlat +datum=WGS84")
408
+ p4s <- CRS("EPSG:4326")
428
409
  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
410
 
430
411
 
@@ -484,45 +465,47 @@
484
465
  ###################################################
485
466
  ### code chunk number 110: die.Rnw:1614-1615
486
467
  ###################################################
487
- writeOGR(Fires[,c("gml_id", "FireDate", "Area_HA")], dsn="fires.kml", layer="fires", driver="KML", overwrite_layer=TRUE)
488
-
468
+ #writeOGR(Fires[,c("gml_id", "FireDate", "Area_HA")], dsn="fires.kml", layer="fires", driver="KML", overwrite_layer=TRUE)
469
+ st_write(st_as_sf(Fires[,c("gml_id", "FireDate", "Area_HA")]), dsn="fires.kml", layer="fires", driver="KML", delete_layer=TRUE)
489
470
 
490
471
  ###################################################
491
472
  ### code chunk number 111: die.Rnw:1673-1679
492
473
  ###################################################
493
- library(maptools)
474
+ #library(maptools)
494
475
  grd <- as(meuse.grid, "SpatialPolygons")
495
476
  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)
477
+ grd.union <- st_union(st_as_sf(grd))
478
+ ll <- CRS("EPSG:4326")
479
+ grd.union.ll <- st_transform(grd.union, st_crs(ll))
499
480
 
500
481
 
501
482
  ###################################################
502
483
  ### code chunk number 112: die.Rnw:1703-1707
503
484
  ###################################################
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))
485
+ library(stars)
486
+ llGRD <- st_as_stars(st_bbox(grd.union.ll), n=600*600)
487
+ llSGDF <- as(st_crop(llGRD, grd.union.ll), "Spatial")
507
488
  llSPix <- as(llSGDF, "SpatialPixelsDataFrame")
508
489
 
509
490
 
510
491
  ###################################################
511
492
  ### code chunk number 113: die.Rnw:1722-1724
512
493
  ###################################################
513
- meuse_ll <- spTransform(meuse, CRS("+proj=longlat +datum=WGS84"))
494
+ meuse_ll <- as(st_transform(st_as_sf(meuse), st_crs(CRS("EPSG:4326"))), "Spatial")
514
495
  llSPix$pred <- gstat::idw(log(zinc)~1, meuse_ll, llSPix)$var1.pred
515
496
 
516
497
 
517
498
  ###################################################
518
499
  ### code chunk number 114: die.Rnw:1744-1749
519
500
  ###################################################
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
501
  image(llSPix, "pred", col=bpy.colors(20))
523
- dev.off()
524
- kmlOverlay(llGRD, "zinc_IDW.kml", "zinc_IDW.png")
525
502
 
526
503
  ###################################################
527
504
  source("die_snow.R", echo=TRUE)
528
505
 
506
+ (sI <- sessionInfo()) # check: no sp?
507
+
508
+ "rgdal" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
509
+ "rgeos" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
510
+ "maptools" %in% c(names(sI$otherPkgs), names(sI$loadedOnly))
511
+