2

I have KML file from downloaded https://data.gov.sg/dataset/electoral-boundary_2020. As you can see on the website, there are 31 sub-polygons and each of them has a label. I have loaded the KML on R.

I have a separate .csv data file with postal codes, latitude and longitude. I would like to extract the name of the polygons from the KML file to map to the .csv to give 4 columns:

postal_code latitude longitude Name

Where Name is the name the polygon if the postal_code fall within the polygon. How do I do this using R?

Here is a sample data

library(rgdal)
library(spatialEco)
library(dplyr)

postal_code=c(117606, 640215, 600220, 628389, 118560, 267314, 129580, 268458)
latitude=c(1.283146, 1.348578, 1.340993, 1.303928, 1.277510, 1.321747, 1.315483, 1.322019)
longitude=c(103.7739, 103.7121, 103.7374, 103.6586, 103.7893, 103.7999, 103.7650, 103.7910)
df=cbind.data.frame(postal_code, latitude, longitude)

coordinates(df)<-~longitude+latitude

poly.sg=readOGR('Downloads/doc.kml', "ELD2020")
proj4string(poly.sg)
poly.sg=spTransform(poly.sg, CRS('+init=epsg:4326'))

I have tried to follow the steps in Join spatial point data to polygons in R

poly.sg2<-as(poly.sg, "SpatialPolygonsDataFrame")
pts.poly<-point.in.poly(df, poly.sg2)

I received this error:

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : st_crs(x) == st_crs(y) is not TRUE

5
  • Have you seen this article gis.stackexchange.com/questions/137621/…?
    – Taras
    Oct 23, 2020 at 6:51
  • @Taras Thank you for sharing the link. I have attempted to follow the steps in the link but encountered an error.
    – tatami
    Oct 23, 2020 at 8:12
  • What kind of error?
    – MarcM
    Oct 23, 2020 at 8:50
  • @MarcM it is in my question
    – tatami
    Oct 23, 2020 at 13:25
  • Without seeing what is being imported from readOGR I can only speculate but, you do not need to coerce to an sp class as, it is already there. Meaning, the polg.sg2 object is not necessary. Besides, the point.in.poly function handles both sp and sf objects. Dec 25, 2021 at 19:28

1 Answer 1

2

You are not setting the CRS to the df object; if you print it you'll notice

class       : SpatialPointsDataFrame 
features    : 8 
extent      : 103.6586, 103.7999, 1.27751, 1.348578  (xmin, xmax, ymin, ymax)
crs         : NA 

That crs is NA; so you have to set it:

library(raster) # for its Spatial*DataFrame print method
library(sp)
#...
coordinates(df)<-~longitude+latitude
sp::proj4string(df) = sp::CRS('+init=epsg:4326')

Then, when you print the df object again, you'll see:

class       : SpatialPointsDataFrame 
features    : 8 
extent      : 103.6586, 103.7999, 1.27751, 1.348578  (xmin, xmax, ymin, ymax)
crs         : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +tow

Now your intersection should work

3
  • Is crs() in crs(df) from the package raster?
    – tatami
    Oct 25, 2020 at 16:30
  • yes, raster::crs(); I've edited the code; also, the object print in the code is from raster package
    – Elio Diaz
    Oct 25, 2020 at 17:03
  • spatialEco imports sp and raster so both proj4strings and crs methods should be available
    – Elio Diaz
    Oct 25, 2020 at 18:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.