.. _chshgeorst:

=============================
Cheat sheet on Geocoordinates
=============================


.. only:: html

    **Links:** :download:`notebook <chsh_geo.ipynb>`, :downloadlink:`html <chsh_geo2html.html>`, :download:`PDF <chsh_geo.pdf>`, :download:`python <chsh_geo.py>`, :downloadlink:`slides <chsh_geo.slides.html>`, :githublink:`GitHub|_doc/notebooks/cheat_sheets/chsh_geo.ipynb|*`


Cheat sheet on geocoding.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





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 <https://pypi.python.org/pypi/pyshp>`__,
`fiona <https://pypi.python.org/pypi/fiona>`__\ …

.. code:: ipython3

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




.. parsed-literal::
    ['.\\Centerline.dbf',
     '.\\Centerline.prj',
     '.\\Centerline.shp',
     '.\\Centerline.shx']



.. code:: ipython3

    import shapefile
    shp = shapefile.Reader("Centerline.shp")

.. code:: ipython3

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

.. code:: ipython3

    shapes[0].__dict__




.. parsed-literal::
    {'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]}



.. code:: ipython3

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




.. parsed-literal::
    {'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 <https://dragons8mycat.wordpress.com/2014/04/29/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.

.. code:: ipython3

    %load_ext pyensae

.. code:: ipython3

    %head Centerline.prj






.. raw:: html

    <pre>
    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]]
    </pre>



This means to be converted into another coordinates system: *ESPG:XXXX*.
That’s what the website `Prj2EPSG <http://prj2epsg.org/search>`__ does.
The module `pyproj <https://pypi.python.org/pypi/pyproj/>`__ does the
conversion into longitude / latitude.

.. code:: ipython3

    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)




.. parsed-literal::
    (-109.78898903640629, 48.55518509715834)