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:
If I plot it with tricontourf, it looks like (errors outlined with red color):
When I plot it with tricontourf and default values it looks like this:
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:
- 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.
- 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
