edzer / sp

Classes and methods for spatial data
http://edzer.github.io/sp/
127 stars 27 forks source link

gridded function problem #62

Open jamessong opened 4 years ago

jamessong commented 4 years ago

Sth. is wrong with the gridded function, is there a way to put the tolerance in it? i tried to used the points2grid function, but I had to divide the value by 1000

gridded(spPnts) <- TRUE Error in if (max(object@grid.index) > .NumberOfCells(object@grid)) return("grid.index max value too large") : missing value where TRUE/FALSE needed In addition: Warning messages: 1: In points2grid(points, tolerance, round) : grid has empty column/rows in dimension 1 2: In points2grid(points, tolerance, round) : grid topology may be corrupt in dimension 1 3: In points2grid(points, tolerance, round) : grid has empty column/rows in dimension 2 4: In points2grid(points, tolerance, round) : grid topology may be corrupt in dimension 2 5: In getGridIndex(coordinates(points), grid) : NAs introduced by coercion to integer range

edzer commented 4 years ago

We cannot do anything here if you don't provide a reproducible example.

jamessong commented 4 years ago
        x        y

1 2.724091 35.25607 2 2.991394 35.27221 3 3.258881 35.28751 4 3.526489 35.30197 5 3.794220 35.31559 6 4.062073 35.32837 7 4.330048 35.34031 8 4.598145 35.35141 9 4.866333 35.36167 10 5.134613 35.37109 11 5.402985 35.37967 12 5.671448 35.38741 13 5.939941 35.39429 14 6.208527 35.40034 15 6.477142 35.40556 16 6.745819 35.40991 17 7.014526 35.41343 18 7.283264 35.41611 19 7.552032 35.41793 20 7.820801 35.41892 21 8.089600 35.41906 22 8.358368 35.41836 23 8.627136 35.41681 24 8.895874 35.41441 25 9.164612 35.41119 26 9.433289 35.40710 27 9.701935 35.40218 28 9.970520 35.39640 29 10.239075 35.38979 30 10.507538 35.38234 31 10.775940 35.37404 32 11.044250 35.36489 33 11.312469 35.35492 34 11.580597 35.34410 35 11.848602 35.33243 36 12.116516 35.31994 37 12.384277 35.30660 38 12.651947 35.29242 39 12.919464 35.27740 40 13.186829 35.26156 41 13.454041 35.24486 42 13.721100 35.22735 43 13.988007 35.20899 44 14.254700 35.18979 45 14.521240 35.16978 46 14.787567 35.14893 47 15.053711 35.12725 48 15.319641 35.10474 49 15.585358 35.08142 50 15.850861 35.05725 51 16.116150 35.03228 52 16.381165 35.00645 53 16.645966 34.97984 54 16.910492 34.95239 55 17.174774 34.92413 56 17.438812 34.89506 57 17.702545 34.86515 58 17.966003 34.83445 59 18.229187 34.80293 60 18.492065 34.77060 61 18.754669 34.73747 62 19.016937 34.70353 63 19.278900 34.66878 64 19.540558 34.63323 65 19.801880 34.59689 66 20.062866 34.55973 67 20.323486 34.52180 68 20.583771 34.48306 69 20.843719 34.44353 70 21.103302 34.40323 71 2.703735 35.47443 72 2.972076 35.49064 73 3.240570 35.50602 74 3.509216 35.52054 75 3.777985 35.53423 76 4.046875 35.54707 77 4.315887 35.55907 78 4.584991 35.57024 79 4.854218 35.58054 80 5.123535 35.59000 81 5.392944 35.59863 82 5.662415 35.60640 83 5.931976 35.61334 84 6.201599 35.61942 85 6.471252 35.62465 86 6.740967 35.62903 87 7.010712 35.63256 88 7.280518 35.63526 89 7.550293 35.63710 90 7.820129 35.63808 91 8.089935 35.63822 92 8.359772 35.63751 93 8.629578 35.63596 94 8.899353 35.63356 95 9.169128 35.63030 96 9.438843 35.62619 97 9.708527 35.62125 98 9.978149 35.61546 99 10.247742 35.60881 100 10.517242 35.60132 101 10.786652 35.59298 102 11.056000 35.58379 103 11.325256 35.57376 104 11.594421 35.56289 105 11.863464 35.55117 106 12.132416 35.53860 107 12.401215 35.52519 108 12.669891 35.51095 109 12.938446 35.49586 110 13.206848 35.47992 111 13.475067 35.46315 112 13.743164 35.44555 113 14.011047 35.42710 114 14.278778 35.40781 115 14.546326 35.38770 116 14.813690 35.36675 117 15.080841 35.34496 118 15.347778 35.32234 119 15.614502 35.29889 120 15.880981 35.27462 121 16.147278 35.24950 122 16.413300 35.22357 123 16.679077 35.19681 124 16.944611 35.16923 125 17.209900 35.14082 126 17.474884 35.11161 127 17.739624 35.08157 128 18.004059 35.05070 129 18.268219 35.01903 130 18.532074 34.98655 131 18.795624 34.95324 132 19.058868 34.91914 133 19.321808 34.88422 134 19.584412 34.84850 135 19.846680 34.81196 136 20.108612 34.77464 137 20.370209 34.73650 138 20.631439 34.69758 139 20.892303 34.65787 140 21.152832 34.61734

The above is a sptial point (convert it to) then do a transformation "+proj=lcc +lon_0=8 +lat_0=54 +lat_1=37 +lat_2=71 +a=6370000.0 +b=6370000.0"

spPntsWRF <- spTransform(spPnts,CRS(projWRF))

gridded(spPntsWRF) <- TRUE

edzer commented 4 years ago

That is an (incomplete) sketch of what you did, but not a reproducible example.

jamessong commented 4 years ago

Here is the Rdata file for your reference On Wednesday, November 6, 2019, 02:20:56 PM EST, Edzer Pebesma notifications@github.com wrote:

That is an (incomplete) sketch of what you did, but not a reproducible example.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

mdsumner commented 4 years ago

Try points2grid with a tolerance, inferring regular grids from points is bad enough but having to guess at the projection and transform and play with tolerance is really bad. I hope you can restore the grid in a better way (I'd prefer to define the grid abstractly, then identify the cell each point belongs to than do this kind of heuristic - but might get you out of a bind).

x <- read.table(text = "1 2.724091 35.25607
2 2.991394 35.27221
3 3.258881 35.28751
4 3.526489 35.30197
5 3.794220 35.31559
6 4.062073 35.32837
7 4.330048 35.34031
8 4.598145 35.35141
9 4.866333 35.36167
10 5.134613 35.37109
11 5.402985 35.37967
12 5.671448 35.38741
13 5.939941 35.39429
14 6.208527 35.40034
15 6.477142 35.40556
16 6.745819 35.40991
17 7.014526 35.41343
18 7.283264 35.41611
19 7.552032 35.41793
20 7.820801 35.41892
21 8.089600 35.41906
22 8.358368 35.41836
23 8.627136 35.41681
24 8.895874 35.41441
25 9.164612 35.41119
26 9.433289 35.40710
27 9.701935 35.40218
28 9.970520 35.39640
29 10.239075 35.38979
30 10.507538 35.38234
31 10.775940 35.37404
32 11.044250 35.36489
33 11.312469 35.35492
34 11.580597 35.34410
35 11.848602 35.33243
36 12.116516 35.31994
37 12.384277 35.30660
38 12.651947 35.29242
39 12.919464 35.27740
40 13.186829 35.26156
41 13.454041 35.24486
42 13.721100 35.22735
43 13.988007 35.20899
44 14.254700 35.18979
45 14.521240 35.16978
46 14.787567 35.14893
47 15.053711 35.12725
48 15.319641 35.10474
49 15.585358 35.08142
50 15.850861 35.05725
51 16.116150 35.03228
52 16.381165 35.00645
53 16.645966 34.97984
54 16.910492 34.95239
55 17.174774 34.92413
56 17.438812 34.89506
57 17.702545 34.86515
58 17.966003 34.83445
59 18.229187 34.80293
60 18.492065 34.77060
61 18.754669 34.73747
62 19.016937 34.70353
63 19.278900 34.66878
64 19.540558 34.63323
65 19.801880 34.59689
66 20.062866 34.55973
67 20.323486 34.52180
68 20.583771 34.48306
69 20.843719 34.44353
70 21.103302 34.40323
71 2.703735 35.47443
72 2.972076 35.49064
73 3.240570 35.50602
74 3.509216 35.52054
75 3.777985 35.53423
76 4.046875 35.54707
77 4.315887 35.55907
78 4.584991 35.57024
79 4.854218 35.58054
80 5.123535 35.59000
81 5.392944 35.59863
82 5.662415 35.60640
83 5.931976 35.61334
84 6.201599 35.61942
85 6.471252 35.62465
86 6.740967 35.62903
87 7.010712 35.63256
88 7.280518 35.63526
89 7.550293 35.63710
90 7.820129 35.63808
91 8.089935 35.63822
92 8.359772 35.63751
93 8.629578 35.63596
94 8.899353 35.63356
95 9.169128 35.63030
96 9.438843 35.62619
97 9.708527 35.62125
98 9.978149 35.61546
99 10.247742 35.60881
100 10.517242 35.60132
101 10.786652 35.59298
102 11.056000 35.58379
103 11.325256 35.57376
104 11.594421 35.56289
105 11.863464 35.55117
106 12.132416 35.53860
107 12.401215 35.52519
108 12.669891 35.51095
109 12.938446 35.49586
110 13.206848 35.47992
111 13.475067 35.46315
112 13.743164 35.44555
113 14.011047 35.42710
114 14.278778 35.40781
115 14.546326 35.38770
116 14.813690 35.36675
117 15.080841 35.34496
118 15.347778 35.32234
119 15.614502 35.29889
120 15.880981 35.27462
121 16.147278 35.24950
122 16.413300 35.22357
123 16.679077 35.19681
124 16.944611 35.16923
125 17.209900 35.14082
126 17.474884 35.11161
127 17.739624 35.08157
128 18.004059 35.05070
129 18.268219 35.01903
130 18.532074 34.98655
131 18.795624 34.95324
132 19.058868 34.91914
133 19.321808 34.88422
134 19.584412 34.84850
135 19.846680 34.81196
136 20.108612 34.77464
137 20.370209 34.73650
138 20.631439 34.69758
139 20.892303 34.65787
140 21.152832 34.61734
")

spPnts <- sp::SpatialPoints(as.matrix(x[c(2,3)]), proj4string = sp::CRS("+proj=longlat +a=6370000.0 +b=6370000.0"))

projWRF <- "+proj=lcc +lon_0=8 +lat_0=54 +lat_1=37 +lat_2=71 +a=6370000.0 +b=6370000.0"
spPntsWRF <- sp::spTransform(spPnts, sp::CRS(projWRF))

sp::points2grid(spPntsWRF, tolerance = 0.0005)
#> Warning in sp::points2grid(spPntsWRF, tolerance = 5e-04): grid has empty
#> column/rows in dimension 1
#>                           V2          V3
#> cellcentre.offset -482727.26 -2012728.26
#> cellsize            24543.02    24542.63
#> cells.dim              70.00        2.00

Created on 2019-11-07 by the reprex package (v0.3.0)

jamessong commented 4 years ago

it is strange that if you divide the spPntWRF by 1000, you can use a 0.002 tolerance (which is fine, 2 meter accuracy). I am so baffled by this ... a bug? spPntsWRF@coords<-spPntsWRF@coords/1000

sp::points2grid(spPntsWRF, tolerance = 0.0015)                   coords.x1   coords.x2cellcentre.offset -482.72851 -2012.72777cellsize            24.54289    24.54381cells.dim           70.00000    66.00000

On Thursday, November 7, 2019, 02:51:00 PM EST, Michael Sumner <notifications@github.com> wrote:  

Try points2grid with a tolerance, inferring regular grids from points is bad enough but having to guess at the projection and transform and play with tolerance is really bad. I hope you can restore the grid in a better way (I'd prefer to define the grid abstractly, then identify the cell each point belongs to than do this kind of heuristic - but might get you out of a bind). x <- read.table(text = "1 2.724091 35.25607 2 2.991394 35.27221 3 3.258881 35.28751 4 3.526489 35.30197 5 3.794220 35.31559 6 4.062073 35.32837 7 4.330048 35.34031 8 4.598145 35.35141 9 4.866333 35.36167 10 5.134613 35.37109 11 5.402985 35.37967 12 5.671448 35.38741 13 5.939941 35.39429 14 6.208527 35.40034 15 6.477142 35.40556 16 6.745819 35.40991 17 7.014526 35.41343 18 7.283264 35.41611 19 7.552032 35.41793 20 7.820801 35.41892 21 8.089600 35.41906 22 8.358368 35.41836 23 8.627136 35.41681 24 8.895874 35.41441 25 9.164612 35.41119 26 9.433289 35.40710 27 9.701935 35.40218 28 9.970520 35.39640 29 10.239075 35.38979 30 10.507538 35.38234 31 10.775940 35.37404 32 11.044250 35.36489 33 11.312469 35.35492 34 11.580597 35.34410 35 11.848602 35.33243 36 12.116516 35.31994 37 12.384277 35.30660 38 12.651947 35.29242 39 12.919464 35.27740 40 13.186829 35.26156 41 13.454041 35.24486 42 13.721100 35.22735 43 13.988007 35.20899 44 14.254700 35.18979 45 14.521240 35.16978 46 14.787567 35.14893 47 15.053711 35.12725 48 15.319641 35.10474 49 15.585358 35.08142 50 15.850861 35.05725 51 16.116150 35.03228 52 16.381165 35.00645 53 16.645966 34.97984 54 16.910492 34.95239 55 17.174774 34.92413 56 17.438812 34.89506 57 17.702545 34.86515 58 17.966003 34.83445 59 18.229187 34.80293 60 18.492065 34.77060 61 18.754669 34.73747 62 19.016937 34.70353 63 19.278900 34.66878 64 19.540558 34.63323 65 19.801880 34.59689 66 20.062866 34.55973 67 20.323486 34.52180 68 20.583771 34.48306 69 20.843719 34.44353 70 21.103302 34.40323 71 2.703735 35.47443 72 2.972076 35.49064 73 3.240570 35.50602 74 3.509216 35.52054 75 3.777985 35.53423 76 4.046875 35.54707 77 4.315887 35.55907 78 4.584991 35.57024 79 4.854218 35.58054 80 5.123535 35.59000 81 5.392944 35.59863 82 5.662415 35.60640 83 5.931976 35.61334 84 6.201599 35.61942 85 6.471252 35.62465 86 6.740967 35.62903 87 7.010712 35.63256 88 7.280518 35.63526 89 7.550293 35.63710 90 7.820129 35.63808 91 8.089935 35.63822 92 8.359772 35.63751 93 8.629578 35.63596 94 8.899353 35.63356 95 9.169128 35.63030 96 9.438843 35.62619 97 9.708527 35.62125 98 9.978149 35.61546 99 10.247742 35.60881 100 10.517242 35.60132 101 10.786652 35.59298 102 11.056000 35.58379 103 11.325256 35.57376 104 11.594421 35.56289 105 11.863464 35.55117 106 12.132416 35.53860 107 12.401215 35.52519 108 12.669891 35.51095 109 12.938446 35.49586 110 13.206848 35.47992 111 13.475067 35.46315 112 13.743164 35.44555 113 14.011047 35.42710 114 14.278778 35.40781 115 14.546326 35.38770 116 14.813690 35.36675 117 15.080841 35.34496 118 15.347778 35.32234 119 15.614502 35.29889 120 15.880981 35.27462 121 16.147278 35.24950 122 16.413300 35.22357 123 16.679077 35.19681 124 16.944611 35.16923 125 17.209900 35.14082 126 17.474884 35.11161 127 17.739624 35.08157 128 18.004059 35.05070 129 18.268219 35.01903 130 18.532074 34.98655 131 18.795624 34.95324 132 19.058868 34.91914 133 19.321808 34.88422 134 19.584412 34.84850 135 19.846680 34.81196 136 20.108612 34.77464 137 20.370209 34.73650 138 20.631439 34.69758 139 20.892303 34.65787 140 21.152832 34.61734 ")

spPnts <- sp::SpatialPoints(as.matrix(x[c(2,3)]), proj4string = sp::CRS("+proj=longlat +a=6370000.0 +b=6370000.0"))

projWRF <- "+proj=lcc +lon_0=8 +lat_0=54 +lat_1=37 +lat_2=71 +a=6370000.0 +b=6370000.0" spPntsWRF <- sp::spTransform(spPnts, sp::CRS(projWRF))

sp::points2grid(spPntsWRF, tolerance = 0.0005)

> Warning in sp::points2grid(spPntsWRF, tolerance = 5e-04): grid has empty

> column/rows in dimension 1

> V2 V3

> cellcentre.offset -482727.26 -2012728.26

> cellsize 24543.02 24542.63

> cells.dim 70.00 2.00

Created on 2019-11-07 by the reprex package (v0.3.0)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

mdsumner commented 4 years ago

Please don't hack S4 slots like that, the bounding box is not updated and so the object is now invalid (amongst other issues).

I don't think there's any bug here, as mentioned above deriving grids from raw points and transform hackery is a desperate workaround at best. If you want help with specific details I suggest to post a reprex with access to the complete set of points and more details on their background on stackoverflow or gis.stackexchange.

jamessong commented 4 years ago

I thought i sent the file, but here it is .

edzer.latlon.Rdata.zip

The above i as actually not using the spPntsWRF@coords<-spPntsWRF@coords/1000 i converted them to the required points for points2grid, why divided by 1000 is fine?

jamessong commented 4 years ago

or maybe the gridded function should have a tolerance in it?

edzer commented 4 years ago

sp is in maintenance mode, gridded does not have a tolerance argument and we are not going to change this. As @mdsumner pointed out, points2grid has a tolerance argument.

Active development is done in sf, and there is a reason why functions like this were not ported to sf: it's too messy.

jamessong commented 4 years ago

glad to hear about that. Thanks.  Now the only question left is that, why when I divide by 1000 i can use 0.0015 tolerance while without dividing by 1000, i can't even use 10 as tolerance? On Thursday, November 7, 2019, 03:40:38 PM EST, Edzer Pebesma notifications@github.com wrote:

sp is in maintenance mode, gridded does not have a tolerance argument and we are not going to change this. As @mdsumner pointed out, points2grid has a tolerance argument.

Active development is done in sf, and there is a reason why functions like this were not ported to sf: it's too messy.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.