search

dieghernan / tidyterra

tidyverse and ggplot2 methods for terra spatial objects
https://dieghernan.github.io/tidyterra/
Other
168 stars 7 forks source link

Warning in geom_spatvector_text #114

Closed rhgof closed 10 months ago

rhgof commented 1 year ago

Am continually getting the warning from geom_spatvector_text

Warning message:
In st_point_on_surface.sfc(sf::st_zm(x)) :
  st_point_on_surface may not give correct results for longitude/latitude data

The data is very simple and for the same data I do not get the error when plotting via geom_spatvector

This is the minimal code that will produce the warning

library(tidyverse)
library(terra)
#> terra 1.7.39
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
library(tidyterra)
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:stats':
#> 
#>     filter

data("world.cities",package = "maps")
citiesOfInterest <- c("Sydney","Wollongong","Newcastle")
cities <- world.cities |>
  filter(country.etc == "Australia") |>
  filter(name %in% citiesOfInterest)

cityVect <- vect(cities,geom=c("long","lat"),crs = "epsg:4326")

ggplot() +
  geom_spatvector(data = cityVect,color="grey40") +
  geom_spatvector_text(data = cityVect,color="grey10",
                       mapping = aes(label=values(cityVect)$name))
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

Created on 2023-09-22 with reprex v2.0.2 The geom_spatvector_text call throws the warning If I comment it out - no warning.

Here is the data

> cities
              name country.etc     pop    lat   long capital
1        Byron Bay   Australia    6812 -28.66 153.61       0
2    Coffs Harbour   Australia   63525 -30.30 153.12       0
3 Forster-Tuncurry   Australia   17718 -32.19 152.52       0
4        Newcastle   Australia  500085 -32.92 151.75       0
5   Port Macquarie   Australia   50421 -31.44 152.91       0
6           Sydney   Australia 4444513 -33.87 151.21       0
7       Wollongong   Australia  261883 -34.42 150.87       0
> cityVect
 class       : SpatVector 
 geometry    : points 
 dimensions  : 7, 4  (geometries, attributes)
 extent      : 150.87, 153.61, -34.42, -28.66  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :             name country.etc   pop capital
 type        :            <chr>       <chr> <int>   <int>
 values      :        Byron Bay   Australia  6812       0
                  Coffs Harbour   Australia 63525       0
               Forster-Tuncurry   Australia 17718       0
dieghernan commented 1 year ago

Yes, this warning comes from sf package. It happens on geom_spatvector_text because that function uses centroids for placing labels, and the warning just highlights that centroids on EPSG:4326 (lon/lat) are not accurate (see explanation on the sf package).

See the (very) simple code of geom_spatvector_text, that is a simple wrapper of geom_sf_text:

https://github.com/dieghernan/tidyterra/blob/eaa55f933344afcd783ba0ce20e93761c33500c1/R/geom_spatvector.R#L123-L143

If you want to avoid the warning because it may be annoying your best option is to project your SpatVector to a suitable CRS, for example:

...

cityVect <- vect(cities,geom=c("long","lat"),crs = "epsg:4326")
# Project to webmercator
cityVect <- project(cityVect, "epsg:3857")

On finding the best CRS for your data I suggest the package crsuggest, that gives you potential CRS candidates depending of the location you are working with.

rhgof commented 1 year ago

Well Australia Datum GDA2020 EPSG:7844 fails with the same issue. Looks like I will have to live with it 

library(tidyverse)
library(terra)
#> terra 1.7.39
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
library(tidyterra)
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:stats':
#> 
#>     filter

data("world.cities",package = "maps")
citiesOfInterest <- c("Sydney","Wollongong","Newcastle")
cities <- world.cities |>
  filter(country.etc == "Australia") |>
  filter(name %in% citiesOfInterest)

cityVect <- vect(cities,geom=c("long","lat"),crs = "EPSG:7844")

ggplot() +
  geom_spatvector(data = cityVect,color="grey40") +
  geom_spatvector_text(data = cityVect,color="grey10",
                       mapping = aes(label=values(cityVect)$name))
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

Created on 2023-09-23 with reprex v2.0.2

dieghernan commented 10 months ago

Seems that ESPG:7844 is non-projected

sf::st_crs("EPSG:7844") |> sf::st_is_longlat()
#> TRUE

If it is only the warning that annoyes you, maybe you want to project your data (introducing crsuggest).

library(raster)
#> Loading required package: sp
library(tidyverse)
library(terra)
#> terra 1.7.55
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
library(tidyterra)
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:raster':
#> 
#>     select
#> The following object is masked from 'package:stats':
#> 
#>     filter

data("world.cities",package = "maps")
citiesOfInterest <- c("Sydney","Wollongong","Newcastle")
cities <- world.cities |>
  filter(country.etc == "Australia") |>
  filter(name %in% citiesOfInterest)

cityVect <- vect(cities,geom=c("long","lat"),crs = "EPSG:7844")

# Maybe project?

crs <- crsuggest::suggest_crs(cityVect, type = "projected")

crs
#> # A tibble: 10 × 6
#>    crs_code crs_name                        crs_type crs_gcs crs_units crs_proj4
#>    <chr>    <chr>                           <chr>      <dbl> <chr>     <chr>    
#>  1 8058     GDA2020 / NSW Lambert           project…    7844 m         +proj=lc…
#>  2 3308     GDA94 / NSW Lambert             project…    4283 m         +proj=lc…
#>  3 7856     GDA2020 / MGA zone 56           project…    7844 m         +proj=ut…
#>  4 28356    GDA94 / MGA zone 56             project…    4283 m         +proj=ut…
#>  5 20256    AGD66 / AMG zone 56             project…    4202 m         +proj=ut…
#>  6 7845     GDA2020 / GA LCC                project…    7844 m         +proj=lc…
#>  7 4462     WGS 84 / Australian Centre for… project…    4326 m         +proj=lc…
#>  8 3577     GDA94 / Australian Albers       project…    4283 m         +proj=ae…
#>  9 3112     GDA94 / Geoscience Australia L… project…    4283 m         +proj=lc…
#> 10 32756    WGS 84 / UTM zone 56S           project…    4326 m         +proj=ut…

ggplot() +
  geom_spatvector(data = cityVect,color="grey40") +
  geom_spatvector_text(data = cityVect,color="grey10",
                       mapping = aes(label=name)) +
  coord_sf(crs=as.numeric(crs$crs_code[1]))

Created on 2023-11-15 with reprex v2.0.2

Not much to do in the tidyterra side, closing it now.