1

Given a geographical coordinate in U.S., how to find out if it is in Urban or Rural areas?

I have about 10000 geographical coordinates all in the U.S., and I want to use Python + basemap to find out if a point is urban or rural.

I'm not sure which library or shape file to use.

I'll need a function like this:

def is_urban(coordinate):
  # use the shapefile
  urban = False
  return urban
Ash
  • 3,428
  • 1
  • 34
  • 44
  • You could for sure use over() methods in R out of the sp package, but I realize that isn't what you asked for (python). The R method would be to find the urban rural classifications, which are likely in the tiger shapefiles from the US census. SpatialPoints your coordinates -> over() using sp. – Neal Barsch Feb 26 '18 at 03:57
  • 1
    Edit to above since I forgot rgdal in R: SpatialPoints your coordinates -> load shapefile using readOGR from rgdal, over() using sp. The tiger shapefiles would still probably be your best bet in python. Hope this helps a little! – Neal Barsch Feb 26 '18 at 04:03
  • 1
    do you have a dataset that delineates urban vs rural areas? – Paul H Feb 28 '18 at 19:52
  • This is the urban areas shape file, so if a point is not inside these, it's rural. https://www.census.gov/geo/maps-data/data/cbf/cbf_ua.html – Ash Feb 28 '18 at 23:47

1 Answers1

5
import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface

pt = (-97.759615,30.258773) # an x,y tuple
shp = shapefile.Reader('/home/af/Downloads/cb_2016_us_ua10_500k/cb_2016_us_ua10_500k.shp') #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()

def is_urban(pt):
    result = False
    for i in range(len(all_shapes)):
        boundary = all_shapes[i] # get a boundary polygon
        #name = all_records[i][3] + ', ' + all_records[i][4] # get the second field of the corresponding record
        if Point(pt).within(shape(boundary)): # make a point and see if it's in the polygon
            result = True
    return result

result = is_urban(pt)

I ended up using shapely and shapefile downloaded from https://www.census.gov/geo/maps-data/data/cbf/cbf_ua.html, which has urban areas of the U.S., so if a point is not within any of these areas, it's rural.

I tested it, and it works to my expectation.

Ash
  • 3,428
  • 1
  • 34
  • 44