[R-sig-Geo] Help with spTransform() function and final plot colors

[Michael Sumner]

Hi, your longlat coordinates are in the western hemisphere (-119, 39) but
EPSG 27700 is designed for (and appropriate for some scenarios) somewhere
in the UK. You need a local projection suitable for the area you are in.
Look for an EPSG for your region, or maybe Albers conic or Lambert
Azimuthal local authority (look on a paper map if you have one). Or define
your own with projstrings (or wkt, it's a bit complicated atm to advise
simple construction). There's some discussion here:
https://geocompr.robinlovelace.net/reproj-geo-data.html#which-crs-to-use

Some will say use UTM, but you should not unless it's a particular local
one used commonly in your work and the region the data is from.

HTH, Mike





On Mon., 18 May 2020, 19:46 Poling, William via R-sig-Geo, <
r-sig-geo using r-project.org> wrote:

> #RStudio Version Version 1.2.1335
> sessionInfo()
> # R version 4.0.0 Patched (2020-05-03 r78349)
> #Platform: x86_64-w64-mingw32/x64 (64-bit) #Running under: Windows 10 x64
> (build 17763)
>
> Hello. I am running my data through a routine I found that finds clusters
> of data points based on distance rule.
>
> https://gis.stackexchange.com/questions/64392/finding-clusters-of-points-based-distance-rule-using-r
>
> 1. I get this error when I get to this point in the routine, see complete
> routine below?
> xy <- spTransform(xy, CRS("+init=epsg:27700 +datum=WGS84")) non finite
> transformation detected:
>
> I tried converting columns from numeric to Integer but did not help.
>
> However, I continue to run the rest of the routine and the final sequence,
> the plot itself, seems to work Can this error be corrected somehow, despite
> the fact that it continues to work, just curious what it is I guess"
>
> 2. However in the plot at the end of the routine the color key appears
> with the colors but the clusters in the plot are black, see plot call at
> the end of the routine below?
>
> Thank you for any advice.
>
> WHP
>
> ERROR:
> non finite transformation detected:
>       coords.x1 coords.x2
>  [1,] -119.7339  39.53939 Inf Inf
>  [2,] -119.7665  39.39630 Inf Inf
>  [3,] -119.7794  39.28768 Inf Inf
>  [4,] -121.0234  39.20503 Inf Inf
>  [5,] -122.0047  47.19262 Inf Inf
>  [6,] -122.0135  47.18883 Inf Inf
>  [7,] -122.0379  47.52190 Inf Inf
>  [8,] -122.0578  47.60975 Inf Inf
>  [9,] -122.1330  47.13669 Inf Inf
> [10,] -122.1509  47.55962 Inf Inf
> [11,] -122.1706  47.15546 Inf Inf
> [12,] -122.1846  47.23485 Inf Inf
> [13,] -122.1846  48.15307 Inf Inf
> [14,] -122.1851  47.44870 Inf Inf
> [15,] -122.1954  47.68485 Inf Inf
> [16,] -122.1990  47.51610 Inf Inf
> [17,] -122.2014  47.44772 Inf Inf
> [18,] -122.2025  47.69815 Inf Inf
> [19,] -122.2037  47.67190 Inf Inf
> [20,] -122.2090  47.40378 Inf Inf
> [21,] -122.2108  47.25336 Inf Inf
> [22,] -122.2291  47.63880 Inf Inf
> [23,] -122.2419  47.76870 Inf Inf
> [24,] -122.2722  48.04803 Inf Inf
> [25,] -122.2732  47.87700 Inf Inf
> [26,] -122.2804  47.77620 Inf Inf
> [27,] -122.2839  47.82103 Inf Inf
> [28,] -122.2890  47.86899 Inf Inf
> [29,] -122.2993  47.67306 Inf Inf
> [30,] -122.3180  47.38217 Inf Inf
> [31,] -122.3270  47.40378 Inf Inf
> [32,] -122.3474  47.43884 Inf Inf
> [33,] -122.3484  47.53083 Inf Inf
> [34,] -122.3581  47.27678 Inf Inf
> [35,] -122.3618  47.76735 Inf Inf
> [36,] -122.3700  47.56567 Inf Inf
> [37,] -122.3908  47.54938 Inf Inf
> [38,] -122.4128  47.64622 Inf Inf
> [39,] -122.4293  47.17660 Inf Inf
> [40,] -122.4621  47.44505 Inf Inf
> [41,] -122.4904  47.27460 Inf Inf
> [42,] -122.5515  46.93979 Inf Inf
> [43,] -122.7348  42.37320 Inf Inf
> [44,] -122.7827  47.31059 Inf Inf
> [45,] -122.7987  47.23475 Inf Inf
> [46,] -122.8385  42.35119 Inf Inf
> [47,] -122.8537  42.34495 Inf Inf
> [48,] -122.8904  42.39555 Inf Inf
> [49,] -122.8927  42.33022 Inf Inf
> [50,] -122.9451  47.37574 Inf Inf
> [51,] -122.9594  42.30376 Inf Inf
> [52,] -123.0641  47.16428 Inf Inf
> [53,] -123.3413  42.44117 Inf Inf
> Error in spTransform(xSP, CRSobj, ...) :
>   failure in points
> 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
> In addition: Warning message:
> In spTransform(xSP, CRSobj, ...) : 53 projected point(s) not finite
>
> Here is the data:
> dput(sample)
> structure(list(ID = 1:53, Longitude = c(-119.733899, -119.766493,
> -119.779416, -121.0234, -122.004736, -122.013456, -122.0379, -122.0578,
> -122.132971, -122.150901, -122.170608, -122.18462, -122.184639,
> -122.185079, -122.195398, -122.198994, -122.201356, -122.202507,
> -122.20371, -122.209047, -122.210797, -122.229095, -122.2419, -122.27216,
> -122.273164, -122.280355, -122.28389, -122.289043, -122.299261,
> -122.318009, -122.326987, -122.347382, -122.34844, -122.358115,
> -122.361839, -122.37003, -122.390815, -122.41282, -122.429323, -122.462136,
> -122.490417, -122.551483, -122.734847, -122.782736, -122.798669,
> -122.838498, -122.853683, -122.8904, -122.89271, -122.94511, -122.959407,
> -123.064087, -123.341346), Latitude = c(39.53939, 39.396298, 39.287681,
> 39.205028, 47.192616, 47.188833, 47.5219, 47.609748, 47.13669, 47.559616,
> 47.155455, 47.234849, 48.15307, 47.448697, 47.684854, 47.516104, 47.447723,
> 47.698146, 47.6719, 47.403778, 47.253364, 47.638795, 47.768701, 48.048027,
> 47.876997, 47.776205, 47.821029, 47.868987, 47.673056, 47.382165,
> 47.403785, 47.438836, 47.530831, 47.276776, 47.76735, 47.565667, 47.549377,
> 47.646222, 47.176596, 47.445053, 47.274599, 46.939789, 42.373195,
> 47.310595, 47.234748, 42.351189, 42.344953, 42.395547, 42.33022, 47.375736,
> 42.303755, 47.164278, 42.441172)), class = "data.frame", row.names = c(NA,
> -53L))
>
> Here is the routine:
> require(sp)
> require(rgdal)
> #25*1609.34 = 40233.5
> dis <- 40233.5 #Distance threshold 25miles converted to meters
>
> x <- tmp1b[,c(6)]
> head(x)
> y <- tmp1b[,c(5)]
> head(y)
> str(y)
>
> xy <- SpatialPointsDataFrame(matrix(c(x,y), ncol=2),
> data.frame(ID=seq(1:length(x))),
>                              proj4string=CRS("+proj=longlat +ellps=WGS84
> +datum=WGS84"))
>
> str(xy)
>
> xy <- spTransform(xy, CRS("+init=epsg:27700 +datum=WGS84"))#Error
> occurring here?
>
> chc <- hclust(dist(data.frame(rownames=rownames(xy using data),
> x=coordinates(xy)[,1], #What is xy using data???? see str(xy)
>                               y=coordinates(xy)[,2])), method="complete")
> str(chc)
>
> view(xy using data)#Look like 3 clusters
>
> # Distance with a 25 Mile threshold
> chc.dis <- cutree(chc, h=dis)
> str(chc.dis)
> view(chc.dis)
>
> # Join results to meuse sp points
> xy using data <- data.frame(xy using data, Clust=chc.dis)
> str(xy using data)
> view(xy using data)
>
> # Plot results
> plot(xy, col=factor(xy using data$Clust), pch=19)
> box(col="black")
> title(main="Clustering")
> legend("topleft", legend=paste("Cluster", 1:3,sep=""),#Change from 4 in
> demo to 3 in my data
>        col=palette()[1:3], pch=rep(19,3), bg="white")
>
>
> Proprietary
>
> NOTICE TO RECIPIENT OF INFORMATION:\ This e-mail may con...{{dropped:16}}
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]