Cheat sheet on Geocoordinates#

Links: notebook, html, PDF, python, slides, GitHub

Cheat sheet on geocoding.

from jyquickhelper import add_notebook_menu
add_notebook_menu()

Read shapefiles#

Shapefiles usually comes in three file:

  • <name>.shp: the shapes

  • <name>.shx

  • <name>.dbf: the metadata

Module you can use to read the files: pyshp, fiona…

from pyquickhelper.filehelper import unzip_files
unzip_files("geo_data.zip", ".")

['.\Centerline.dbf',
 '.\Centerline.prj',
 '.\Centerline.shp',
 '.\Centerline.shx']
import shapefile
shp = shapefile.Reader("Centerline.shp")

records = shp.records()
shapes = shp.shapes()
fields = shp.fields

shapes[0].__dict__

{'shapeType': 3,
 'points': [(1312867.164001003, 260498.0600337088),
  (1313098.4133638293, 260489.8426001817),
  (1313098.288386017, 260482.84261640906)],
 'parts': [0],
 'bbox': [1312867.164001003, 260482.84261640906, 1313098.4133638293, 260498.0600337088]}

{k[0]:v for k,v in zip(fields[1:], records[0])}

{'StreetName': 'NE 118th CT',
 'FromStreet': '132nd PL NE',
 'ToStreet': 'E of 132nd PL NE',
 'StreetWidt': 34.0,
 'ConstYear': 1980.0,
 'ConstBy': '',
 'd_Classifi': 4.0,
 'd_Ownershi': 'RED',
 'd_Status': 'ACT',
 'd_DataSour': 'EXT',
 'd_InsideCi': 'YES',
 'FromLeft': 13200,
 'ToLeft': 13398,
 'FromRight': 13201,
 'ToRight': 13399,
 'LeftCity': 'RED',
 'RightCity': 'RED',
 'LeftZip': '98052',
 'RightZip': '98052',
 'Side': 'B',
 'Exclude': '',
 'Location': '',
 'lz_evn': 'G',
 'lz_odd': 'G',
 'la_evn': 'G',
 'la_odd': 'G',
 'StreetNa_1': 'NE 118 CT',
 'FireArteri': 0,
 'surface': '1',
 'Migrated': 'YES',
 'd_ConstBy': '',
 'RedCenterl': 7,
 'MaxSpeedLi': 25,
 'CreatedBy': '',
 'DateCreate': b'',
 'ModifiedBy': 'ESP-NFEARS',
 'DateModifi': datetime.date(2016, 7, 21),
 'StAlias1': '',
 'StAlias2': '',
 'StAlias3': '',
 'AssetID': 'Cntl7',
 'LucityAuto': 0,
 'd_AssetTyp': 'Centerline',
 'CalcLength': 238.39641915,
 'Shape_len': 238.396419151}

Convert coordinates into longitude / latitude#

Source : GIS Tips – How to Find the EPSG Code of your Shapefile

The file <name>.prj contains informations about the coordinates used to encode the shapefiles.

%load_ext pyensae

%head Centerline.prj

PROJCS["NAD_1983_HARN_StatePlane_Washington_North_FIPS_4601_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-120.8333333333333],PARAMETER["Standard_Parallel_1",47.5],PARAMETER["Standard_Parallel_2",48.73333333333333],PARAMETER["Latitude_Of_Origin",47.0],UNIT["Foot_US",0.3048006096012192]],VERTCS["NGVD_1929",VDATUM["National_Geodetic_Vertical_Datum_1929"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Foot_US",0.3048006096012192]]

This means to be converted into another coordinates system: ESPG:XXXX. That’s what the website Prj2EPSG does. The module pyproj does the conversion into longitude / latitude.

from pyproj import Proj, transform
p1 = Proj(init='epsg:2926')  # returned by Prj2ESPG
p2 = Proj(init='epsg:4326')  # longitude / latidue
x0, y0 = 1312487.225652799, 231314.96255882084
transform(p1,p2,x0,y0)

(-109.78898903640629, 48.55518509715834)