4

I have a shapefile of US county boundaries which I got from the NHGIS. I can load the file using the rgdal package and see information about it:

library(rgdal)
shp_1860 <- readOGR("shp", "US_county_1860", stringsAsFactors = F)
summary(shp_1860)

# Output:
Object of class SpatialPolygonsDataFrame
Coordinates:
       min     max
x -2356112 2258225
y -1337373 1565782
Is projected: TRUE 
proj4string :
[+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0]

You'll notice that the projection is a little unusual. But what I really don't understand is the latitude and longitude: what do those values mean? Why aren't they latitudes and longitudes in the vicinity of the United States?

I can fortify this shapefile for use in ggplot2, and plot it without a problem.

library(ggplot2)
shp_1860_df <- fortify(shp_1860, region="GISJOIN")
p <- ggplot() +
  geom_map(map = shp_1860_df, data= shp_1860_df,
           color = "black", fill = "white", size = 0.1,
           aes(map_id = id)) +
  expand_limits(x = shp_1860_df$long, y = shp_1860_df$lat) +
  coord_equal()
print(p)

I get a map just fine with that code. If I replaced coord_equal() with the more logical coord_map(), R and my computer always crashes because memory use skyrockets. I presume this is because of the weird coordinates.

I thought I would re-project the shapefile in order to plot it.

shp_1860_proj <- spTransform(shp_1860, CRS("+proj=longlat +datum=WGS84"))
summary(shp_1860_proj)

# Output
Object of class SpatialPolygonsDataFrame
Coordinates:
        min       max
x -124.7625 -66.94993
y   24.5211  49.38436
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]

Notice that the latitude and longitude now have values that make sense. I can plot this new shapefile using base R without a problem, e.g., plot(shp_1860_proj). But when I try to fortify it, I get an error:

shp_1860_proj_df <- fortify(shp_1860_proj, region="GISJOIN")
Loading required package: rgeos
rgeos version: 0.3-1, (SVN revision 413M)
 GEOS runtime version: 3.3.3-CAPI-1.7.4 
 Polygon checking: TRUE 

Error in createPolygonsComment(p) : 
  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 6

What is the right way to go about re-projecting and then fortifying this shapefile? Alternatively, can I re-project it after I've fortified it?

  • 1
    You may find partly working solution in this thread. – SS_Rebelious Oct 24 '13 at 20:41
  • Also try out this solution. Seems that there is something wrong with the ring of the one of your polygons. – SS_Rebelious Oct 24 '13 at 20:48
  • They aren't "long" and "lat", despite what fortify calls them. They are X, and Y (or less meaningfully "Easting" and "Northing"), in a metres-based coordinate system as defined by the PROJ.4 string. – mdsumner Aug 30 '16 at 13:42
  • The rgeos error comes from fortify's attempts to union the polygons on "GISJOIN", which is more interesting. Do you need to dissolve on that attribute? What if you just leave out the region argument to fortify? Also, can you post your R version as well? rgeos 0-3-1 is very old now. – mdsumner Aug 30 '16 at 13:45

Your Answer

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

Browse other questions tagged or ask your own question.