What projections in r will fatten a city map?

Question

I'm trying to make some city maps, and I'd like to distort them a bit so that the densest parts (manhattan in this case) appear to be a little wider on the final map. I'm having a hard time finding good projection examples at the city-level.

I'd just like a projection that makes a better-looking map of New York City, but that might be too subjective, so let's say I want to fatten Manhattan. Reproducible code of what I've tried below:

library(ggplot2)
library(maptools)

shpny.tf <- tempfile()

download.file(
    "http://www2.census.gov/geo/tiger/TIGER2010/COUNTY/2010/tl_2010_36_county10.zip" ,
    shpny.tf ,
    mode = 'wb'
)

shpny.uz <- unzip( shpny.tf , exdir = tempdir() )

ny.shp <- readShapePoly( shpny.uz[ grep( 'shp$' , shpny.uz ) ] )

# limit the shapefile to only the five boroughs
nyc.shp <- subset( ny.shp , as.numeric( as.character( COUNTYFP10 ) ) %in% c( 5 , 47 , 61 , 81 , 85 ) )

# prepare for plotting
nyc <- fortify( nyc.shp )

# make a plot
p <- qplot( data = nyc , x = long , y = lat )

# default plot
p

# these two seem like good candidates for fattening,
# but i'm not sure how to make them take up more space on the right and left side
p + coord_map("ortho", orientation=c(41, -74, -25))
p + coord_map("azequalarea",orientation=c(41,-74,-22))

# this one's okay but still not ideal
p + coord_map("lagrange")

# this one is named for new york yet i can't get it to do anything useful
p + coord_map("newyorker",r=10)

ideally, the result would look closer to the proportions in this new york subway map. i understand that this much of a distortion isn't possible with only projection, i'm just curious if any projections could make the map closer to that ideal shape. thanks!

Answer 1

I downloaded the shapefile and stored it since constantly downloading a 1M file has questionable efficacy. Since I did that with curl on the command line, I was already in a terminal, so I pre-extracted the five boroughs with ogr2ogr:

$ ogr2ogr -f "ESRI Shapefile" \ 
          -where 'COUNTYFP10 IN ("005", "047", "061", "081", "085")' \ 
          boroughs.shp tl_2010_36_county10.shp

Next, I read the new shapefile in with readOGR and setup some necessary libraries:

library(rgdal)
library(maptools)
library(ggplot2)

boroughs.shp <- readOGR("tl_2010_36_county10/", "boroughs")

Now, we separate New York and the other four boroughs:

nyc <- boroughs.shp[boroughs.shp$NAME10 == "New York",]
other_four <- boroughs.shp[boroughs.shp$NAME10 != "New York",]

Read up on ?elide as it lets you both scale polygons and move them around and you'll need to do more than I'm showing in this example. Here, I scale nyc by 1.1 (arbitrarily) but keep the 1:1 scale for the others. I'm taking a bit of a shortcut here and doing both since elide changes the CRS for the SpatialPolygonsDataFrames.

nyc <- elide(nyc, scale=1.1)
other_four <- elide(other_four, scale=1)

Next, we re-combine them:

boroughs.shp <- spRbind(nyc, other_four)

and do the necessary work for ggplot:

boroughs <- fortify(boroughs.shp, region="NAME10" )

I prefer geom_map and this give you contiguous polygons vs your point plot:

gg <- ggplot(data=boroughs, aes(x=long, y=lat))
gg <- gg + geom_map(map=boroughs, aes(map_id=id, group=group), fill="white", color="black")
gg <- gg + coord_map()
gg

That, obviously, is not the end-state you want, but you've now got the tools you need to build a proper cartogram (though I am personally not a fan of distorted polygon maps).

You can read up more on making cartograms to get the final effect, but it'll take some more effort. Real cartograms are the way to go here, IMO, if you want to do the fattening/scaling to show density.

EDIT OR, perhaps just settle for globular?

gg_globular <- gg + coord_map("globular", orientation=c(41,-74,-22)) + ggtitle("globular")

Answer 2

Still in the vein of cartograms, you could combine R with the cross-platform Java program, ScapeToad. Clearly this isn't as desirable as a pure-R solution, but it's an option all the same.

For example:

# Create a field that will be used to distort polygons. 
# In our case, the second polygon corresponds to New York County (i.e. Manhattan)
ny.shp$size <- c(1, 1.5, 1, 1, 1)

Now duck over to ScapeGoat and generate the cartogram, using the size field as the cartogram attribute, and export as .shp.

This produces the following:

Clearly this approach is less appropriate if you need accurate grid lines.

Answer 3

Pseudocode for a solution in whatever language you are most comfortable writing in. This approach would allow you to apply different levels of distortion to different parts of the map whilst preserving the topologies of common edges.

I suspect messing around with different projections unlikely to produce any noticeable distortion over such a small geographic area.

Anyway, this is how I would approach it.

1. Dump points of each polygon, along with poly_id and index position.
2. Pick a reference point. For example, centroid of your densest polygon
3. Calculate azimuth and distance between each poly point (step 1) and your reference point (step 3).
4. Create a translation function with whatever parameters you desire. It could be simply a function of distance keeping azimuth constant or distance and sin of azimuth etc. if you want to distort more in a particular plane.
5. Calculate the new x,y of each poly point
6. Remake the polygons.

I can already think of a great data animation showing the shapes morphing over time. When I get some time, I'll code up the solution and share. HTH

Answer 4

http://kartoweb.itc.nl/geometrics/coordinate%20transformations/coordtrans.html

Refer the theory. :-(((

My recommendation would be to do all your analysis using your current CRS.

I would then warp your results by vectorizing every layer including the rasters (eg dump pixel centroids). Next I'd apply a translation function to distort features of every vector object in the same manner. Then I would your re-rasterize your Krigged point estimates.

My tool of choice for this would be postgis using ST_translate. I'd then pull the data back into R to create new raster surfaces for visualization purposes.

Answer 5

# this is clearly insane
p +
    geom_map(map=nyc, aes(map_id=id, group=group), fill="white", color="black") +
    coord_map( "globular" , orientation = c( 40.55 , -74.11 , -24 ) )