Use a projected coordinate system in the states.json example file

[Nowosad]

Hey @jbaileyh, I would suggest to use a projected coordinate system in the states.json example file. A comparison between a geographic and projected coordinate system is attached below. You can read more about projections at https://source.opennews.org/articles/choosing-right-map-projection/.

library(geogrid)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.4, proj.4 4.9.3
input_file = system.file("extdata", "states.json", package = "geogrid")
us = st_read(input_file)
#> Reading layer `OGRGeoJSON' from data source `/home/jn/R/x86_64-redhat-linux-gnu-library/3.4/geogrid/extdata/states.json' using driver `GeoJSON'
#> Simple feature collection with 49 features and 5 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -124.7332 ymin: 24.5447 xmax: -66.9499 ymax: 49.38436
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
plot(st_geometry(us))

us_2163 = st_transform(us, 2163)
plot(st_geometry(us_2163))

Created on 2018-05-01 by the reprex package (v0.2.0).

[jbaileyh]

Hey @Nowosad this is an interesting suggestion. Originally I planned on keeping clear of projections (leaving it as +proj=longlat +datum=WGS84 +no_defs) and enabling the user to define whatever is most appropriate for their use-case.

Other than aesthetics (visually I much prefer regionally appropriate projections), do you have any reasons in mind for why to do this? Given your background it would be great to hear your perspective / best practice recommendations.

Side note: I need to check but this change may require a modification to the use of spDists depending on the projection of the underlying data.

[Nowosad]

My thinking:

1. Using projected data shows people a best practice.
2. The output looks better (see below).
3. I've got spDist warnings on the original data (a lot of them, I just left a few below as an example). Projected version does not give warnings.

My side question: Why it is important to have two separate functions - calculate_grid and assign_polygons? I would assume that users just want to transform their dataset (e.g. an sp object -> the same object as hexagons). Have you though about having just one function, for example transform_to_grid()?

        library(geogrid)
        library(sf)
        library(tmap)
        input_file = system.file("extdata", "states.json", package = "geogrid")
        us = st_read(input_file)        
        us_2163 = st_transform(us, 2163)

        hex_cells1 = calculate_grid(us, grid_type = "hexagonal", seed = 25, learning_rate = 0.03)
        us_hex1 = assign_polygons(us, hex_cells1)
#> Warning in sp::spDists(original_points, new_points, longlat = FALSE):
#> spDists: argument longlat conflicts with CRS(x); using the value FALSE

#> Warning in sp::spDists(original_points, new_points, longlat = FALSE):
#> spDists: argument longlat conflicts with CRS(x); using the value FALSE

        hex_cells2 = calculate_grid(us_2163, grid_type = "hexagonal", seed = 25, learning_rate = 0.03)
        us_hex2 = assign_polygons(us_2163, hex_cells2)

        tm1 = tm_shape(us_hex1) + tm_polygons()
        tm2 = tm_shape(us_hex2) + tm_polygons()

        tmap_arrange(tm1, tm2)

Created on 2018-05-01 by the reprex package (v0.2.0).

[sassalley]

@Nowosad submitting to CRAN but having issues with valgrind which is new to me - wondering if you might have insight?

https://www.stats.ox.ac.uk/pub/bdr/memtests/valgrind/geogrid/tests/testthat.Rout

[Nowosad]

Sadly no. It looks like something related to rcpp... Maybe @hafen can help?