2

I'm trying to interpolate data within boundaries and plot contour lines (polygons) based on Latitude, longitude, value data from csv files.

Results I want print as geojson.

I'm stuck on the basic contour plot from csv. I really appreciate help here.

This is what I got in this moment but can't make it work.

import numpy as np
import matplotlib.pyplot as plt


data = np.genfromtxt('temp.csv', delimiter=',')

x = data[:1]
y = data[:2]
[x,y] = meshgrid(x,y);
z = data[:3];

plt.contour(x,y,z)
plt.show()

The CSV looks like this

2015-10-28T09:25:12.800Z,51.56325,17.529043333,4.6484375,19.8046875
2015-10-28T09:25:12.900Z,51.56325,17.529041667,4.4921875,19.6875
2015-10-28T09:25:13.000Z,51.563248333,17.529041667,4.453125,19.9609375
2015-10-28T09:25:13.200Z,51.563245,17.529041667,4.140625,19.765625
2015-10-28T09:25:13.300Z,51.563243333,17.529041667,4.609375,19.6484375
2015-10-28T09:25:13.500Z,51.563241667,17.529041667,4.609375,19.53125
2015-10-28T09:25:13.600Z,51.56324,17.529041667,4.4921875,19.375
2015-10-28T09:25:13.700Z,51.563238333,17.529041667,4.4140625,19.765625
2015-10-28T09:25:13.800Z,51.563236667,17.529041667,4.453125,20.234375
2015-10-28T09:25:13.900Z,51.563235,17.529041667,4.3359375,19.84375
2015-10-28T09:25:14.000Z,51.563233333,17.529041667,4.53125,19.453125
2015-10-28T09:25:14.100Z,51.563231667,17.529041667,4.53125,19.8046875
2015-10-28T09:25:14.200Z,51.563228333,17.529041667,4.1796875,19.4921875
2015-10-28T09:25:14.300Z,51.563226667,17.529041667,4.2578125,19.453125
2015-10-28T09:25:14.400Z,51.563225,17.529041667,4.4921875,19.4921875
2015-10-28T09:25:14.500Z,51.563223333,17.529041667,4.375,19.453125
2015-10-28T09:25:14.600Z,51.563221667,17.529041667,4.609375,18.90625
2015-10-28T09:25:14.700Z,51.563218333,17.529041667,4.53125,19.6875
2015-10-28T09:25:14.900Z,51.563215,17.529041667,4.140625,18.75
2015-10-28T09:25:15.000Z,51.563213333,17.529041667,4.453125,18.828125

Column 1 - Latitude Column 2 - Longitude Column 3 - Value

For contour lines I need also limits - for example (4.1,4.3,4.6)

seek
  • 1,065
  • 16
  • 33
  • I guess there is a few step between your code and having a geojson as output. Did you experience some error message when running your code ? – mgc Jan 19 '16 at 23:11

2 Answers2

9

I guess there is some mistake in your code (according to your data you shouldn't do x = data[:1] but more x = data[..., 1]).

With your of data, the basic steps I will follow to interpolate the z value and fetch an output as a geojson would require at least the shapely module (and here geopandas is used for the convenience).

import numpy as np
import matplotlib.pyplot as plt
from geopandas import GeoDataFrame
from matplotlib.mlab import griddata
from shapely.geometry import Polygon, MultiPolygon

def collec_to_gdf(collec_poly):
    """Transform a `matplotlib.contour.QuadContourSet` to a GeoDataFrame"""
    polygons, colors = [], []
    for i, polygon in enumerate(collec_poly.collections):
        mpoly = []
        for path in polygon.get_paths():
            try:
                path.should_simplify = False
                poly = path.to_polygons()
                # Each polygon should contain an exterior ring + maybe hole(s):
                exterior, holes = [], []
                if len(poly) > 0 and len(poly[0]) > 3:
                    # The first of the list is the exterior ring :
                    exterior = poly[0]
                    # Other(s) are hole(s):
                    if len(poly) > 1:
                        holes = [h for h in poly[1:] if len(h) > 3]
                mpoly.append(Polygon(exterior, holes))
            except:
                print('Warning: Geometry error when making polygon #{}'
                      .format(i))
        if len(mpoly) > 1:
            mpoly = MultiPolygon(mpoly)
            polygons.append(mpoly)
            colors.append(polygon.get_facecolor().tolist()[0])
        elif len(mpoly) == 1:
            polygons.append(mpoly[0])
            colors.append(polygon.get_facecolor().tolist()[0])
    return GeoDataFrame(
        geometry=polygons,
        data={'RGBA': colors},
        crs={'init': 'epsg:4326'})


data = np.genfromtxt('temp.csv', delimiter=',')
x = data[..., 1]
y = data[..., 2]
z = data[..., 3]
xi = np.linspace(x.min(), x.max(), 200)
yi = np.linspace(y.min(), y.max(), 200)
zi = griddata(x, y, z, xi, yi, interp='linear') # You could also take a look to scipy.interpolate.griddata

nb_class = 15 # Set the number of class for contour creation
# The class can also be defined by their limits like [0, 122, 333]
collec_poly = plt.contourf(
    xi, yi, zi, nb_class, vmax=abs(zi).max(), vmin=-abs(zi).max())

gdf = collec_to_gdf(collec_poly)
gdf.to_json()
# Output your collection of features as a GeoJSON:
# '{"type": "FeatureCollection", "features": [{"type": "Feature", "geometry": {"type": "Polygon", "coordinates": [[[51.563214073747474,
# (...)

Edit: You can grab the colors values (as an array of 4 values in range 0-1, representing RGBA values) used by matplotplib by fetching them on each item of the collection with the get_facecolor() method (and then use them to populate one (or 4) column of your GeoDataFrame :

colors = [p.get_facecolor().tolist()[0] for p in collec_poly.collections]
gdf['RGBA'] = colors

Once you have you GeoDataFrame you can easily get it intersected with your boundaries. Use these boundaries to make a Polygon with shapely and compute the intersection with your result:

mask = Polygon([(51,17), (51,18), (52,18), (52,17), (51,17)])
gdf.geometry = gdf.geometry.intersection(mask)

Or read your geojson as a GeoDataFrame:

from shapely.ops import unary_union, polygonize
boundary = GeoDataFrame.from_file('your_geojson')
# Assuming you have a boundary as linestring, transform it to a Polygon:
mask_geom = unary_union([i for i in polygonize(boundary.geometry)])
# Then compute the intersection:
gdf.geometry = gdf.geometry.intersection(mask_geom)
mgc
  • 5,223
  • 1
  • 24
  • 37
  • Thank you, this will help me a lot. I was try to plot the results first but I have an error if I use all of my data about 40k points. Error in griddata (RuntimeError: In initialize: Triangulation is invalid). Maybe I need to sort this data first ? But how? – seek Jan 20 '16 at 09:57
  • I have boundaries as a Polygon (sometimes even polygon with holes) in geojson file. geometry.intersection gives me very strange results. Showing only smallest polygons inside. About RGBA colors I have error: "ValueError: Length of values does not match length of index" at gdf['RGBA'] = colors – seek Jan 20 '16 at 21:54
  • I'm sorry. I doesn't use this unary_union function on my boundaries polygon. Now its looks like it works.. If you can check why this colors doesnt work for me will be great. – seek Jan 20 '16 at 22:01
  • If your boundaries file has multiple features you should try `mask_geom = unary_union(your_boundary_gdf.geometry)` before `gdf.geometry.intersection(mask_geom)`. About the colors i guess it is because some polygons from the matplotlib collection where not fetched as shapely geometry (invalid ones) so it could explain the mismatch of the length. A safer way would have been to fetch them in the `collec_to_gdf` function for each polygon created. – mgc Jan 20 '16 at 22:01
  • @user3476099 i edited the `collec_to_gdf` function so it should outputs a `GeoDataFrame` with a column for RGBA values. – mgc Jan 20 '16 at 22:12
  • Great. Thank you. But I have the same error as earlier ValueError: Length of values does not match length of index – seek Jan 20 '16 at 22:16
  • You're welcome, i hoped this would work.. So..look at the new edit in this function code, i moved the `colors.append(...)` part, it has more chance now to fit the length of the polygons created. (you also can try to `print(len(colors))` and `print(len(polygons))` to see i there is a big difference (and try to guess from where it comes!)) – mgc Jan 20 '16 at 22:22
  • Works like a charm. Once again thank you! I changed RGBA attribute into two - fill & fill-opacity to fit into geojson guidelines. – seek Jan 20 '16 at 22:36
  • I figure it out that if my masking polygon have hole inside this geometry.intersection doesn't use this hole. Do you know how to handle this? – seek Jan 21 '16 at 15:22
  • How would one do this without `geopandas`, using shapely and pandas ++? – M.T Apr 13 '16 at 12:58
  • @M.T You can totally do the same without geopandas (as it is almost not involved in this code, just as a contained for the result). At the end of the function `collec_to_gdf` you just have to replace `return GeoDataFrame(...` by something like `return polygons` if you want a list (it's a python `list` filled by the shapely polygons previously computed) or by something like `return pd.DataFrame({"my_geoms": polygons})` if you want a pandas DataFrame, etc.. – mgc Apr 13 '16 at 18:01
1

Almost I got what I excepted. Based on mgc answer. I change griddata as you suggest to scipy.interpolate.griddata and interpolation method to nearest. Now its works like I want.

Only last thing that I need is to limit this interpolation to polygon (boundaries) also from geoJson.

The other problem is to export to geojson polygons WITH colors as attributes.

import numpy as np
import matplotlib.pyplot as plt
#from matplotlib.mlab import griddata
from scipy.interpolate import griddata

from geopandas import GeoDataFrame
from shapely.geometry import Polygon, MultiPolygon

def collec_to_gdf(collec_poly):
    """Transform a `matplotlib.contour.QuadContourSet` to a GeoDataFrame"""
    polygons = []
    for i, polygon in enumerate(collec_poly.collections):
        mpoly = []
        for path in polygon.get_paths():
            try:
                path.should_simplify = False
                poly = path.to_polygons()
                # Each polygon should contain an exterior ring + maybe hole(s):
                exterior, holes = [], []
                if len(poly) > 0 and len(poly[0]) > 3:
                    # The first of the list is the exterior ring :
                    exterior = poly[0]
                    # Other(s) are hole(s):
                    if len(poly) > 1:
                        holes = [h for h in poly[1:] if len(h) > 3]
                mpoly.append(Polygon(exterior, holes))
            except:
                print('Warning: Geometry error when making polygon #{}'
                      .format(i))
        if len(mpoly) > 1:
            mpoly = MultiPolygon(mpoly)
            polygons.append(mpoly)
        elif len(mpoly) == 1:
            polygons.append(mpoly[0])
    return GeoDataFrame(geometry=polygons, crs={'init': 'epsg:4326'})

data = np.genfromtxt('temp2.csv', delimiter=',')
x = data[..., 1]
y = data[..., 2]
z = data[..., 4]
xi = np.linspace(x.min(), x.max(), num=100)
yi = np.linspace(y.min(), y.max(), num=100)

#zi = griddata(x, y, z, xi, yi, interp='nn') # You could also take a look to scipy.interpolate.griddata
zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='nearest')

nb_class = [5,10,15,20,25,30,35,40,45,50] # Set the number of class for contour creation
# The class can also be defined by their limits like [0, 122, 333]
collec_poly = plt.contourf(
    xi, yi, zi, nb_class, vmax=abs(zi).max(), vmin=-abs(zi).max())

gdf = collec_to_gdf(collec_poly)
#gdf.to_json()
print gdf.to_json()

plt.plot(x,y)

plt.show()
seek
  • 1,065
  • 16
  • 33