1

When I create GeoDataFrame from source data and then show it by pyplot all good, because I didn't do anything with the data. It looks like this:

Image 1

If I plot it with tricontourf, it looks like (errors outlined with red color):

Image 2

When I plot it with tricontourf and default values it looks like this:

Image 3

All source data are stored in .xlsx format(3 different files: Latitudes, Longitudes and Values). Projection of data is WGS84 and I convert it to Lambert Conformal Conic (Only in geopandas). Attached test files to process (download link below). DOWNLOAD TEST FILES HERE

I see only 2 ways to solve this problem:

  1. Crop result with base Latitudes/Longitudes coordinates with creating some polygon and crop it with creating 2 images (result image and mask from polygon), but after some operation I have wrong polygon which delete some negative and positive data.
  2. Do some process to optimize output result with geopandas or matplotlib (I think need to make some operations with geopandas).

Here is my code:

import csv
import io

import geopandas
import numpy
import pandas
from PIL import Image
from cartopy import crs
from matplotlib import pyplot
from pyproj import CRS
from scipy.interpolate import griddata


def read_file(filename):
    result = list()

    with open(filename) as csv_handler:
        csv_reader = csv.reader(csv_handler)

        for line in csv_reader:
            for item in line[0].split(';'):
                item = float(item)

                if item == -999.0:
                    result.append(numpy.nan)
                elif item == -999.7999877929688:
                    result.append(numpy.nan)
                else:
                    result.append(item)

    return result


def load_data():
    longitudes = read_file('lon_NO20_20201019_08_15_06.csv')
    latitudes = read_file('lat_NO20_20201019_08_15_06.csv')
    values = read_file('tmp_NO20_20201019_08_15_06_300.csv')

    return longitudes, latitudes, values


def create_figure_axes():
    figure = pyplot.figure(figsize=(90, 60))
    axes = pyplot.axes(projection=crs.LambertConformal(standard_parallels=(17, 67), central_longitude=80))
    axes.set_xlim([-4250000.000000000, 4750000.000000000])
    axes.set_ylim([-1000000.000000000, 5000000.000000000])
    axes.set_axis_off()

    return figure, axes


def create_geo_data_frame(longitudes: list, latitudes: list, values: list):
    geo_data_frame = geopandas.GeoDataFrame(
        data=pandas.DataFrame(
            {
                'values': values
            }
        ),
        geometry=geopandas.points_from_xy(longitudes, latitudes),
        crs=CRS('EPSG:4326')
    )
    geo_data_frame.to_crs(crs.LambertConformal(standard_parallels=(17, 67), central_longitude=80).proj4_init, inplace=True)

    return geo_data_frame


def show_image(figure):
    with io.BytesIO() as BytesIO:
        figure.savefig(BytesIO, format='png', bbox_inches='tight', pad_inches=0, transparent=True)
        image = Image.open(BytesIO)
        image.show()


def default(longitudes: list, latitudes: list, values: list):
    figure, axes = create_figure_axes()

    geo_data_frame = create_geo_data_frame(longitudes, latitudes, values)

    geo_data_frame.plot(ax=axes, column='values')

    show_image(figure)

    pyplot.close(figure)


def interpolate(longitudes: list, latitudes: list, values: list):
    points = numpy.array(list(zip(longitudes, latitudes)))

    min_longitudes = min(longitudes)
    max_longitudes = max(longitudes)
    min_latitudes = min(latitudes)
    max_latitudes = max(latitudes)

    quality = 12
    longitudes_quality = complex((max_longitudes - min_longitudes) * quality)
    latitudes_quality = complex((max_latitudes - min_latitudes) * quality)
    grid_longitudes, grid_latitudes = numpy.mgrid[min_longitudes:max_longitudes:longitudes_quality, min_latitudes:max_latitudes:latitudes_quality]
    grid_values = griddata(points, numpy.array(values), (grid_longitudes, grid_latitudes), method='cubic')

    new_longitudes = list()
    for element in grid_longitudes:
        for j in element:
            new_longitudes.append(j)

    new_latitudes = list()
    for element in grid_latitudes:
        for j in element:
            new_latitudes.append(j)

    new_values = list()
    for element in grid_values:
        for j in element:
            new_values.append(j)

    return new_longitudes, new_latitudes, new_values


def with_interpolate(longitudes: list, latitudes: list, values: list):
    figure, axes = create_figure_axes()

    longitudes, latitudes, values = interpolate(longitudes, latitudes, values)

    geo_data_frame = create_geo_data_frame(longitudes, latitudes, values)

    geo_data_frame.plot(ax=axes, column='values')

    show_image(figure)

    pyplot.close(figure)


def tricontourf(longitudes: list, latitudes: list, values: list):
    figure, axes = create_figure_axes()

    geo_data_frame = create_geo_data_frame(longitudes, latitudes, values)

    longitudes = list()
    latitudes = list()
    for point in geo_data_frame.geometry:
        longitudes.append(point.x)
        latitudes.append(point.y)
    values = geo_data_frame['values']
    axes.tricontourf(longitudes, latitudes, values)

    show_image(figure)

    pyplot.close(figure)


def main():
    longitudes, latitudes, values = load_data()

    default(longitudes, latitudes, values)
    with_interpolate(longitudes, latitudes, values)
    tricontourf(longitudes, latitudes, values)


if __name__ == '__main__':
    main()

Update

I still cant do process with creating mask, I didnt understand how to set coefficent

I douing this:

longitudes = numpy.array(self.data[satellite][date][time]['longitudes'])
latitudes = numpy.array(self.data[satellite][date][time]['latitudes'])
cond01 = (longitudes - 1) ** 2 + latitudes ** 2 <= 1
cond02 = (longitudes - 0.7) ** 2 + latitudes ** 2 > 0.3
x = longitudes[cond01 & cond02]
y = latitudes[cond01 & cond02]

triangulation = Triangulation(x, y)

triangles = triangulation.triangles

xtri = x[triangles] - numpy.roll(x[triangles], 1, axis=1)
ytri = y[triangles] - numpy.roll(y[triangles], 1, axis=1)
maxi = numpy.max(numpy.sqrt(xtri ** 2 + ytri ** 2), axis=1)

self.data[satellite][date][time]['mask'] = maxi > 0.1

And it didnt work for me well, it almost didnt work and write error message like

ValueError: x and y arrays must have a length of at least 3

1 Answer 1

1

It is a matplotlib problem, not a geospatial problem.

If you look at Contour plot of irregularly spaced data

  1. contour and contourf expect the data to live on a regular 2D grid. Therefore, you need to obtain a regular grid from the irregular points then the contour plot is based on this 2d grid.

  2. As the names tricontour and tricontourf implies, these functions will perform a triangulation internally and the contour plot is based on this triangulation:

Result of the triangulation (triangles in black, points in red):

enter image description here

Then the contour plot is based on this triangulation and not only the points as you wish.

enter image description here

The solution is to play with the parameters of tricontour to try to resolve the problem

mask(ntri,) array-like of bool, optional
Which triangles are masked out.

as in matplotlib contour/contourf of concave non-gridded data

Result for the W (left) side for example

enter image description here

2
  • Thank you @gene, can you help me again, i add some question about setting mask
    – tonysdev
    Apr 14, 2021 at 2:35
  • @tonysdev Please open a new question rather than editing your question with a follow-up question.
    – Aaron
    Apr 14, 2021 at 2:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.