-
Notifications
You must be signed in to change notification settings - Fork 50
/
Copy pathunproj.Rout.save
96 lines (91 loc) · 3.63 KB
/
unproj.Rout.save
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
R version 4.3.3 (2024-02-29) -- "Angel Food Cake"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> # DOI 10.1007/s11004-011-9344-7
> # http://mypage.iu.edu/~srobeson/Pubs/variogram_sphere_mathgeo_2011.pdf
>
> suppressPackageStartupMessages(library(sp))
> library(gstat)
>
> if (require(sp, quietly = TRUE) &&
+ suppressPackageStartupMessages(require(fields, quietly = TRUE)) &&
+ suppressPackageStartupMessages(require(sf, quietly = TRUE))) {
+ data(meuse)
+ coordinates(meuse) = ~x+y
+ proj4string(meuse) = CRS("+init=epsg:28992")
+ ll = "+proj=longlat +ellps=WGS84"
+ # meuse.ll = spTransform(meuse, CRS("+proj=longlat +ellps=WGS84"))
+ meuse.ll = as(st_transform(sf::st_as_sf(meuse), sf::st_crs(ll)), "Spatial")
+ meuse.ll[1:10,]
+ variogram(log(zinc)~1, meuse.ll)
+
+ cloud1 = variogram(log(zinc)~1, meuse, cloud=T, cutoff=6000)
+ cloud2 = variogram(log(zinc)~1, meuse.ll, cloud=T, cutoff=6)
+
+ plot(cloud1$dist/1000, cloud2$dist, xlab="Amersfoort, km", ylab = "Long/lat")
+ abline(0,1)
+
+ data(ozone2)
+ oz = SpatialPointsDataFrame(ozone2$lon.lat,
+ data.frame(t(ozone2$y)),
+ proj4string=CRS("+proj=longlat +ellps=WGS84"))
+ variogram(X870731~1,oz[!is.na(oz$X870731),])
+ utm16 = "+proj=utm +zone=16"
+ # oz.utm = spTransform(oz, utm16)
+ oz.utm = as(sf::st_transform(sf::st_as_sf(oz), utm16) , "Spatial")
+ variogram(X870731~1,oz.utm[!is.na(oz$X870731),])
+
+ # Timothy Hilton, r-sig-geo, Sept 14, 2008:
+
+ foo <-
+ structure(list(z = c(-1.95824831109744, -1.9158901643563, 4.22211761150161,
+ 3.23356929459598, 1.12038389231868, 0.34613850821113, 1.12589932643631,
+ 23.517912251617, 3.0519158690268, 3.20261431141517, -2.10947106854739
+ ), lon = c(-125.29228, -82.1556, -98.524722, -99.948333, -104.691741,
+ -79.420833, -105.100533, -88.291867, -72.171478, -121.556944,
+ -89.34765), lat = c(49.87217, 48.2167, 55.905833, 56.635833,
+ 53.916264, 39.063333, 48.307883, 40.0061, 42.537756, 44.448889,
+ 46.242017)), .Names = c("z", "lon", "lat"), row.names = c(NA,
+ -11L), class = "data.frame")
+
+ coordinates(foo) <- ~lon+lat
+
+ proj4string(foo) <- CRS('+proj=longlat +ellps=WGS84')
+
+ vg.foo <- variogram(z~1, foo, cloud=TRUE, cutoff=1e10)
+
+ cat('==========\nvariogram:\n')
+ print(head(vg.foo))
+
+ cat('==========\nspDistsN1 Distances:\n')
+ print(spDistsN1(coordinates(foo), coordinates(foo)[1,], longlat=TRUE))
+ }
==========
variogram:
dist gamma dir.hor dir.ver id left right
1 3115.1190 8.971063e-04 0 0 var1 2 1
2 1908.0468 1.909846e+01 0 0 var1 3 1
3 1405.6211 1.883757e+01 0 0 var1 3 2
4 1837.5667 1.347749e+01 0 0 var1 4 1
5 1522.8765 1.325847e+01 0 0 var1 4 2
6 119.9386 4.886139e-01 0 0 var1 4 3
==========
spDistsN1 Distances:
[1] 0.0000 3115.1190 1908.0468 1837.5667 1481.6293 3775.6751 1480.4477
[8] 3081.7541 4090.4022 665.9348 2683.1516
Warning message:
In CPL_crs_from_input(x) :
GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
>
> proc.time()
user system elapsed
1.466 1.314 1.282