R Leaflet – How to Read Geo-Spatial Data

gdalleafletpolygonrrgdal

I have some site boundary data as follow:

UNB_BDSB_SiBdry_v08_201028lico_25831.cpq
UNB_BDSB_SiBdry_v08_201028lico_25831.dbf
UNB_BDSB_SiBdry_v08_201028lico_25831.prj
UNB_BDSB_SiBdry_v08_201028lico_25831.sbn
UNB_BDSB_SiBdry_v08_201028lico_25831.sbx
UNB_BDSB_SiBdry_v08_201028lico_25831.shp
UNB_BDSB_SiBdry_v08_201028lico_25831.shp.xml
UNB_BDSB_SiBdry_v08_201028lico_25831.shx

Reading those is R via rgdal package:

myfile = readOGR(dsn = dsn
                 layer = "UNB_BDSB_SiBdry_v08_201028lico_25831")

I'm getting a warning:

OGR data source with driver: ESRI Shapefile 
Source: dsn, layer: "UNB_BDSB_SiBdry_v08_201028lico_25831"
with 1 features
It has 4 fields
Warning message:
In OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS,  :
  Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

And the data looks like:

head(myfile@data)
            NAME                              GISFILE SHAPE_Leng SHAPE_Area
0 Norfolk Boreas UNB_BDSB_SiBdry_v08_201028lico_25831   90513.37  442457610

Now if I use the addpolygons function from leaflet, I can not see any added polygon!

leaflet::leaflet(data = myfile) %>% addTiles() %>%
addPolygons(
            fillColor = "green",
            highlight = highlightOptions(color = "red", weight = 5,
                                         fillOpacity = 0.7,
                                         bringToFront = TRUE))

I'm wondering if there is something wrong with the data or I am doing something wrong!

Also the result using the st_read function from the sf package is as follow :

sf_data = st_read(dsn = "\\\\sv52hcn8934.corp.vattenfall.com/bgf93$/Documents/Map/220420_AH_BdyShp",
+         layer = "UNB_BDSB_SiBdry_v08_201028lico_25831")


summary(sf_data)
     NAME             GISFILE            SHAPE_Leng      SHAPE_Area                 geometry
 Length:1           Length:1           Min.   :90513   Min.   :442457610   POLYGON      :1  
 Class :character   Class :character   1st Qu.:90513   1st Qu.:442457610   epsg:25831   :0  
 Mode  :character   Mode  :character   Median :90513   Median :442457610   +proj=utm ...:0  
                                       Mean   :90513   Mean   :442457610                    
                                       3rd Qu.:90513   3rd Qu.:442457610                    
                                       Max.   :90513   Max.   :442457610

Best Answer

Here's a sample sp-class points data frame in a projected coordinate system like your data is:

> ps
class       : SpatialPointsDataFrame 
features    : 4 
extent      : 1, 4, 1, 4  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
variables   : 0

If I view it in leaflet, its not in the right place...

> leaflet(ps) %>% addTiles() %>% addMarkers()

and leaflet doesn't complain. But if I convert it to an sf spatial data frame and use leaflet:

> leaflet(st_as_sf(ps)) %>% addTiles() %>% addMarkers()
Warning messages:
1: sf layer is not long-lat data 
2: sf layer has inconsistent datum (+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs).
Need '+proj=longlat +datum=WGS84'

...I get those warnings. This is telling me that leaflet needs data in lat-long coordinates, so I can transform them to coordinate system 4326:

> leaflet(st_transform(st_as_sf(ps),4326)) %>% addTiles() %>% addMarkers()

and that gets me the points in the right place with no errors.

I don't know why leaflet doesn't warn with an sp class object but does with the equivalent sf object - its another good argument for dropping sp objects and using sf as much as possible.

Related Question