[GIS] Finding Statistical Area (SA) in which address (Lat\Long) lies via R Point in Polygon

australiacoordinate systemgoogle mapsmapsr

Australian Bureau of Statistics (ABS) provide census information for all Australia. However, ABS has its own division of a state/suburb/town called Statistical Areas Level 1-4 (SA), which is different from the usual government defined suburb-like divisions. ABS also provides ESRI shapefiles for its SA1 divisions.

Problem at hand is that I have number of addresses and I used google maps to find Lat/Longs for them. Now I want to find out in R-Language as in which SA1 division each address belongs to?

I am following the approach mentioned at this Bear-in-Park problem page. Main question that I am struggling with is that I believe Google map's lat/long are in different CRS than ABS ones. I don't know which CRS system Google map uses but some information about ABS shapefile is as follows, 

"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

So how to assign a CRS to Google-Map downloaded Lat/Longs and then how to transform them in same CRS as ABS and then how to determine which address goes to which sa1?

Best Answer

As you can see in this answer: Which CRS to use for Google Maps?, Google Maps uses EPSG 4326. In proj4 format is +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs. So, in R the procedure for any projection is:

library(sp)

wgs84 <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '

pts <- cbind(c(-29.333,-29.444,-29.555), c(148.555,148.666,148.777)) # dummy point

pts_ <- SpatialPoints(pts, proj4string = CRS(wgs84))

pts_repro <- spTransform(pts_,CRS("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"))

But, the output coordinates are almost the same:

ptos_repro
## SpatialPoints:
##      coords.x1 coords.x2
## [1,]   -29.333   148.555
## [2,]   -29.444   148.666
## [3,]   -29.555   148.777
## Coordinate Reference System (CRS) arguments: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +no_defs 

coordinates(pts_repro) == coordinates(pts_)
##      coords.x1 coords.x2
## [1,]      TRUE     FALSE
## [2,]     FALSE      TRUE
## [3,]      TRUE      TRUE

Both ellipsoids are pretty similar.

If this doesn't work, maybe coordinates are in web mercator (if you don't post examples coordinate we can help you more)

Related Question