import os
import numpy as np
import matplotlib.pyplot as plt
import shapefile as shp
from pprint import pprint

"""
Explanation: <figure>
  <IMG SRC="../../logo/logo.png" WIDTH=250 ALIGN="right">
</figure>

IHE Python course, 2017
Reading, writing and working with shape files
T.N.Olsthoorn
2017-03-15
When working with GIS you import and export shape files. A shape file consists of 3 files with the same name but different extension. The one with .shp exteions contains the actual shapes, coordinates and spatial topology, while the one with the .dbf extension contains a table with the user data that is associated with each shape in the shapefile. In it, each shape file has a record of data values that belong to textual headers that indicated what the values mean. The contents of the .dbf file is, in fact a table, with headers, in which each line represents the data pertaining to one shape in the file. The third file, with .shx extension is merely a list of pointers that allow fast searching the shapefiles in GIS. It.s not important for us.
In fact, shape files are the most widely used means to exchange spatial data.
Shape files can be very useful to spatially attribute values to a model, like the conductivities, land use etc. Knowing how they work and how to deal with them can provide a lot of functionality in dealing with spatial data, with or without a GIS.
If you the comman
import shapefile as shp
does not work, you should install pyshp
fire up your cmd window and type
pip install pyshp
The documentation is in the wiki of the pyshp project on github:
https://github.com/GeospatialPython/pyshp/wiki
The more general functionalities are impoted next:
End of explanation
"""


from matplotlib.path import Path as Polygon

"""
Explanation: We need a function that will tell us whether a coordinate pair is inside a shape or not. This function exists in matplotlib.path and is called Path. Path creates a polygon which can then check if points are inside or outside it. Because I like the word Polygon better than Path for this goal, I import the Path class and call it Polygon here.
End of explanation
"""


def inpoly(x, y, pgcoords):
    """Returns bool array [Ny, Nx] telling which grid points are inside polygon

    """
    try:
        isinstance(pgcoords,(list, tuple, np.ndarray))
        len(pgcoords[0])==2
        pgon = Polygon(pgcoords)
    except:
        print('pgcoords must be like [(0, 0), (1, 0), ..] or\n'
            +'an np.ndarray of shape [Np, 2]')
        raise TypeError("Can't create polygon, pgcoords error")

    try:
        x = np.array(x, dtype=float)
        y = np.array(y, dtype=float)
        x.shape == y.shape
    except:
        raise err.ShapeError("x and y not np.ndarrays with same shape.")

    if len(x.shape)==1:
        Y, Y = np.meshgrid(x, y)
    else:
        X = x
        Y = y
    xy = np.vstack((X.ravel(), Y.ravel())).T
    return pgon.contains_points(xy).reshape(X.shape)

"""
Explanation: I define a function inpoly that checks whether the coordinates x,y are inside the shape given by pgcoords. The function returns a logical array (boolean area) of the same shape as the inputs x and y such that points inside the polygon are True in the array and False otherwise.
Notice that inpoly uses the imported Polygon constructor on line 8.
The function does some import checking and guarantees that the shape of the boolean output array is the same as the input array x and y.
The checking is done in lline 27 using the method pgon.contains_points(xy)
This function is very practical, so you may want to use it in other projects as well.
End of explanation
"""


os.listdir('.')

"""
Explanation: CD to the exercises/Mar21 dictory, where we have a shape file stored in the subdirectory Sectie (Section)
End of explanation
"""


# directory (this will be different on your computer)

shpdir =  os.path.join('.', 'Sectie')

fname  = os.listdir(shpdir)[0]

print("The shapefile to work with: ",fname)

"""
Explanation: Define the director where the shape file is and get the file name
End of explanation
"""


rdr = shp.Reader(os.path.join(shpdir, fname))

print("\nAttributes and methods accessible throuhgh this reader:\n")

# learn to read and create comprehensions, this is one
pprint([p for p in dir(rdr) if not p.startswith('_')])

"""
Explanation: Open a reader to read the shapefile:
Then show its attributes, so that we can see what we can do with it.
End of explanation
"""


rdr.fields

"""
Explanation: So we can read the bounding box of all the shapes in the file, its table, its elevaton (if specified), its fiels names etc.
Let's get some useful information about the shapefile. Here, the names of the data fields of each shape (again a comprehension is used to get them in a list)
End of explanation
"""


fldNames = [p[0] for p in rdr.fields]

print(fldNames)

"""
Explanation: Each field not only gives the name of the data but also whether it string )'C' or a number 'N' and its length (first number). The second number is the number of didgets, when the number is floating point.
To get only the field name in a list, we use an other comprehension:
End of explanation
"""


print(fldNames)
for r in rdr.iterRecords():
    print(r)

"""
Explanation: you may use

rdr.iterShapes() to iterate over the shapes in the shapefile
rdr.iterRecords() to iterate over the records in the shapefile
rdr.iterShapeRecords() to iterate over both simultaneously

Print the records in the dbf file, one for each shape:
End of explanation
"""


shprecs = rdr.shapeRecords()

shprecs

"""
Explanation: Read the shapeRecords from the shapefile into a list:
End of explanation
"""


for i, sr in enumerate(shprecs):
    print("Shape number ", i+1)
    print("Bounding box = ",sr.shape.bbox)
    print("Record       = ",sr.record)
    print()

"""
Explanation: Each of these ShapeRecods contain both the shape and its record (data).
Show that we can itereate over the ShapeRecords, get the shape and the record from it and do somthing with it:
End of explanation
"""


[att for att in dir(sr.shape) if not att.startswith('_')]

"""
Explanation: What attributes/methods has a shape?
Simply use a comprehension to show their names:
End of explanation
"""


rdr.bbox

"""
Explanation: Let's now get serious and define a model grid containing all shapes. Simply use the overall bounding box to get the grid extension.
Here is the box:
End of explanation
"""


ix = [0, 2, 2, 0, 0]  # x-indices in rdr.box
iy = [1, 1, 3, 3, 1]  # y-indices

"""
Explanation: Get the indices from the bounding box so that we can plot is in a single line:
End of explanation
"""


plt.plot(np.array([rdr.bbox[i] for i in ix]),
         np.array([rdr.bbox[i] for i in iy]), 'r', lw=3)
plt.show()

"""
Explanation: With these indices in sequence we can get the x and y coordinates of the bbox so arranged that the bbox is plotted (in red with linewidth 3, so that it shows up at the boundary of the plot.
End of explanation
"""


np.array(shprecs[0].shape.points)

"""
Explanation: Now see how we can access the point os a shape, here the first shape in our shape-records list, so shprecs[0].shape.
End of explanation
"""


plt.plot(np.array([rdr.bbox[i] for i in ix]),
         np.array([rdr.bbox[i] for i in iy]), 'r', lw=3)

for sr in shprecs:
    pnts = np.array(sr.shape.points)
    plt.plot(pnts[:,0], pnts[:,1])
plt.show()

"""
Explanation: With all this in place we can plot the bounding box of all shapes and then each shape itself using its points as an array.
End of explanation
"""


# grid line coordinates
x = np.linspace(rdr.bbox[0], rdr.bbox[2], 81)
y = np.linspace(rdr.bbox[1], rdr.bbox[3], 41)

# cell center coordinates
xm = 0.5 * (x[:-1] + x[1:])
ym = 0.5 * (y[:-1] + y[1:])

# generage a full 2D array for both the x and y coordinates
XM, YM = np.meshgrid(xm, ym)

"""
Explanation: The next step is to fill the grid points that fall within each shape with the data value that belongs to the shape. We'll use the horizontal conductivity KH that is contained in each shape's record as can be seen from the field names above.
First generate a grid. A grid is an array with the x coordinate of each cell and one with the y coordinate of each cell.
We'll choose x, y as the coordinates of the cell boundaries (the grid lines) and xm, ym as the coordinates of the cell centers.
Choosing 81 x-grid lines between the horizontal extend 0 and 80 m gives cells of 1 m width.
Choosing 41 y-grid lines between the vertical exten -10 and 0 m gives cells of 0.25 m height.
End of explanation
"""


fig, ax = plt.subplots()

ikh = fldNames.index('KH')
KH = np.zeros_like(XM)      # 2D array to fill with conductivities from shapes

# iterate ove the list of shaperecords. In each loop the next index,
# the next color for the plot and the next records is given. 
for i, clr, sr in zip(range(len(shprecs)), "brgmcy", shprecs):
    pnts = sr.shape.points             # the shape coordinates
    inarray = inpoly(XM, YM, pnts)     # boolean array to match the shape
    KH[inarray] = sr.record[ikh]       # fill in the value
    ax.plot(XM[inarray], YM[inarray], clr+'o', label="shape {}".format(i))
ax.legend(loc='best', fontsize='x-small')
plt.show()

"""
Explanation: Then fill in the data value for each grid and plot the result, once color per shape.
ikh is the index for the field 'KH' (see fields above)
End of explanation
"""


# plt.spy shows where an array is non-zero

plt.spy(KH)
plt.show() 

"""
Explanation: Et voilà !
Now the array can be used as input for a groundwater model for example.
Some ways to show the array
End of explanation
"""


indices = range(0, KH.shape[1], 8)

for i, kh in zip(indices, np.split(KH, indices[1:], axis=1)):
    print("Columns {}:{}".format(i,i+kh.shape[1]))
    for L in kh:
        print(("{:8.4g}" * len(L)).format(*L))
    print()


"""
Explanation: Show the KH array the Matlab way, which is much easier to see. In spyder you may inspect it in the variable window as a kind of spreadsheet.
Because the array is far too wide to print on one line, we chop it in chunk of 8 columns
and print one chunk after the other.
First generate the indices to cut the array in pieces then use np.split( )
Put this together with the indices in a zip, so that the in each for loop we get
the next index and the next array.
For each chunk plot which columns it has.
Then cut the chunk in lines and print each line separately.
Notice that we generate the format dynamically in the print statement using the
length L, i.e. the number of values in the line.
End of explanation
"""


print("The kh values of the different shapes: ", [sr.record[ikh] for sr in shprecs])

plt.matshow(KH)
plt.show()

"""
Explanation: Another way to show an array is to plot it using plt.matshow. Each cell will be colored according to its value. Clearly if the values differ too little, the show up as the same value. This is actually the case here.
Notice that matshow does not use the coordinates, it uses the index of the rows and columns. Because of this, it's plot upside down, just as was the case with plt.spy above.
End of explanation
"""


plt.contourf(XM, YM, KH, 50)
plt.show()

"""
Explanation: You could also contour the shapes, but that has the same problems as with matshow: because of the small differences between the values of different shapes. However, the contourf function uses the coordinates of the grid, so it's plotted correctly on scale and right side up.
End of explanation
"""


