I am trying to geocode a dataset using the given latitude and longitude. I am then going to use these points and join polygon features to them. However, the projection I am using for the points is incorrect. How do I go about this in a better more efficient way?

I begin by uploading the files and setting the projection. The polygon file I have already has a projection so I used that projection for the points.

fema <- readOGR(paste(fema.path, "/2010_DFIRM.shp", sep=""))
#CRS arguments:
 #+proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 
 #+lon_0=-99 +x_0=600000 +y_0=3999999.9998984 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
 #+units=us-ft +no_defs 

data <- read.table(paste(inf.path, "/samp_HT_ANALYSIS", sep = ""), sep=",", fill=T)
coords <- data[,c(22:23)]
data.sp <- SpatialPointsDataFrame(data = data, coords = coords, proj4string =
            CRS("+proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=3999999.9998984 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"))

However, when I do the spatial join using the following code...

  test <- over(data.sp, fema, returnList = T)

No attributes overlap with my points.

I went on to map the points in ArcMap with a base layer map and they aren't mapped at all on the plane. Can someone please explain how I can solve this? Also, general information on projections in R would be great. I thought I could just set the CRS equal of one file equal to the other but I may have to reproject all of my files.

UPDATE: I have rerun my code using suggestions below. However, they are still not plotting correctly. I checked some of the indivudal lat long coordinates in my data and they are within my study region. Here is the updated code.

# Data temp for map
data.use <- read.csv(paste(inf.path, "/Houston_sample_subset.csv", sep = ""))
coords <- as.data.frame(matrix(0,nrow=nrow(data.use),2))
coords[,1] <- data.use$PROP_LATITUDE
coords[,2] <- data.use$PROP_LONGITUDE
data.sp <- SpatialPointsDataFrame(data = data.use, coords = coords, proj4string =
                                    CRS("+proj=longlat +datum=WGS84 +no_defs"))

# Study area
study.area <- readOGR("./FEMA_FLOOD/tl_2017_us_county/tl_2017_us_county.shp")
study.area <- spTransform(study.area, crs(data.sp))
study.area <- subset(study.area, study.area$CBSAFP == 26420)

tiff("./PLOTS/test.tiff", width = 5, height = 5, units = "in", res = 800)
plot(data.sp, add =T)

The resulting plot is Houston CBSA

Only the study region is plotted. I switched to plotting study region vs. fema maps because the fema maps take too long to generate.


2 Answers 2


I believe this is happening due to your SpatialPointsDataFrame being put into a wrong CRS. Your coordinates are in lat long correct? If so they should correspond to the WGS 1984 prj4, which is:

+proj=longlat +datum=WGS84 +no_defs

What you are doing is assigning lat long coordinates to a different projection. Unlike ArcMap, R does not do on the fly conversions. After assigning the proper coordinates to your spatial points you will have to convert one of your object's CRS to the other using the following function:

Here I am converting the fema shapefile's CRS to WGS 1984 (same projection as your data.sp)

fema <- spTransform(fema, crs(data.sp))

After both your spatial objects are in the same CRS, you can retry your Spatial Join.

Let me know if that makes sense!

  • Where are you seeing the CRS for data.sp? Good answer otherwise.
    – ThetaFC
    Dec 3, 2019 at 19:45
  • DPek had provided data.sp's CRS. They are trying to assign lat long CRS to something else. It should be this: data.sp <- SpatialPointsDataFrame(data = data, coords = coords, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
    – Aliyan Haq
    Dec 3, 2019 at 20:20
  • I only see the CRS provided for fema
    – ThetaFC
    Dec 3, 2019 at 20:53
  • @ThetaFC data.sp doesn't have a CRS given it is just a dataset with lat and long. I assigned the CRS when using SpatialPointsDataFrame to be equal to that of fema
    – DPek
    Dec 3, 2019 at 22:01
  • 1
    @DPek then you should do what @AliyanHaq suggests. Assigning a CRS does not change the projection. It only changes the metadata. You need to align one of the two datasets with the other, using the function spTransform. This function actually transforms the underlying coordinates. Also, a geographic coordinate system can be defined using different datums, with WGS84 and NAD83 being the two most widespread in usage. In other words, you can't necessarily just assume it's WGS84
    – ThetaFC
    Dec 3, 2019 at 22:14

I determined the solution. Lat and long were mixed up. Longitude and lattitude is the correct order (xy).

coords[,2] <- data.use$PROP_LATITUDE
coords[,1] <- data.use$PROP_LONGITUDE

is the correct way to set up the data.

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.