2

I have a list of addresses that I have geocoded to get the corresponding longitude and latitude.

> AddressDetail <- geocode(match_address)
> coordinates(AddressDetail) <- c("lon", "lat")

Now I believe I should assign a projection so that R knows how the lon/lat coords are set up

> proj4string(AddressDetail) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

Now I read in the shapefile that uses the NZMG projection

> Shape_MB = readOGR("/Work in progress/MB06_LV2/", "MB06_LV2")  

Now this is where it starts to go wrong. I try transforming the projection of the lon/lat addresses into NZMG projection.

> AddressDetail <- spTransform(AddressDetail_v2, CRS("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs")                                                      
non finite transformation detected:
lon lat  
Error in spTransform(xSP, CRSobj, ...) : failure in points 
In addition: Warning message:
In spTransform(xSP, CRSobj, ...) : 9 projected point(s) not finite

So when I continue on to use the over function to find the corresponding census meshblock index of the polygon that the address is within I get the following error.

> temp <- over(AddressDetail_v2,Shape_MB,MB06)
Error: identicalCRS(x, y) is not TRUE

I'm a bit of a novice with geospatial stuff in R so would appreciate any advice!

  • It seems that you are using NZTM instead of the NZMG definition. See the first answer to this question. (My apologies - I would normally put this as a comment, but I am not allowed to comment yet.) – cengel Sep 10 '13 at 14:44
  • Check the coordinates in your "AddressDetail" data. I suspect that you have a few points, 9 to be exact, that are not defined correctly (e.g., missing "-"). – Jeffrey Evans Jan 13 '14 at 18:46
1

Now that I have significantly more R experience under my belt, it's time to answer this question...

First was to make sure the geocoded address was converted to a spatialpoint cleaner:

AddressDetail <- geocode(match_address)
coordinates(AddressDetail) <- ~lon + lat

Then the relevant projection was assigned to the spatialpoint

wgs84 <- '+proj=longlat +datum=WGS84'
projection(AddressDetail) <- wgs84

From here I wanted to project the geocoded point to the same projection of the meshblocks I want to match the points over.

So read in the meshblock shapefile:

Shape_MB <- readOGR("/Work in progress/MB06_LV2/", layer = "MB06_LV2")

Then simply copy the projection of the meshblock shapefile when transforming the geocoded points:

AddressDetail <- spTransform(AddressDetail, crs(projection(Shape_MB)))

So it really doesn't matter whether the meshblock shapefile is NZMG or something else, just copy the projection and then projection everything else based on that.

Thanks, Alistair

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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