Python and Shapefile: Very large coordinates after importing shapefile
Question: I downloaded a shapefile of Boston and wants to plot it out using the code
below. However it's giving me an error `ValueError: lat_0 must be between
-90.000000 and 90.000000 degrees `

Turns out `coords` has values `(33869.92130000144, 777617.2998000011,
330800.31099999696, 959741.1853)` Why is it so large?

[Boston shapefile is obtained here](http://www.mass.gov/anf/research-and-
tech/it-serv-and-support/application-serv/office-of-geographic-information-
massgis/datalayers/zipcodes.html)

**Code**

    
    
    # Import Boston shapefile
    shapefilename = 'ZIPCODES_NT_POLY'
    shp = fiona.open(shapefilename + '.shp')
    coords = shp.bounds
    shp.close()
    
    w, h = coords[2] - coords[0], coords[3] - coords[1]
    extra = 0.01
    
    m = Basemap(
        projection='tmerc', ellps='WGS84',
        lon_0 = np.mean([coords[0], coords[2]]), 
        lat_0 = np.mean([coords[1], coords[3]]), 
        llcrnrlon = coords[0] - extra * w,
        llcrnrlat = coords[1] - extra * h,
        urcrnrlon = coords[2] + extra * w,
        urcrnrlat = coords[3] + extra * h,
        resolution = 'i', suppress_ticks = True)
    

**Error**

    
    
    ValueError: lat_0 must be between -90.000000 and 90.000000 degrees
    


Answer: You need to reproject from the native projection of Eastings and Northings to
another coordinate reference system. If you want degrees of Latitude and
Longitude, which is usually
[WGS84](https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84) or
[EPSG:4326](http://epsg.io/4326). Here's how to reproject it:

    
    
    import fiona
    import pyproj
    from functools import partial
    from shapely.geometry import box
    from shapely.ops import transform
    
    shp = fiona.open('ZIPCODES_NT_POLY.shp', 'r')
    p_in = pyproj.Proj(shp.crs)
    bound_box = box(*shp.bounds)
    shp.close()
    
    p_out = pyproj.Proj({'init': 'EPSG:4326'})  # aka WGS84
    project = partial(pyproj.transform, p_in, p_out)
    
    bound_box_wgs84 = transform(project, bound_box)
    
    print('native box: ' + str(bound_box))
    print('WGS84 box: ' + str(bound_box_wgs84))
    

> native box: POLYGON ((330800.310999997 777617.2998000011, 330800.310999997
> 959741.1853, 33869.92130000144 959741.1853, 33869.92130000144
> 777617.2998000011, 330800.310999997 777617.2998000011))
>
> WGS84 box: POLYGON ((-69.93980848410942 41.23787282321487, -69.899038698261
> 42.87724537285449, -73.53324195423785 42.8704709990465, -73.48147096070339
> 41.2312695091688, -69.93980848410942 41.23787282321487))

Otherwise, most of the parameters that are required by Basemap are in
`shp.crs` (take a look).

