0

My ultimate goal is to create geoJSON from GRIB2 format using Python in order to display it at the Leaflet-open-street map.

I have the following data extracted from GRIBs:

Multi-dimensional array of x's (lon)

[[  0.     1.25   2.5  ... 356.25 357.5  358.75]  
 ...  
 [  0.     1.25   2.5  ... 356.25 357.5  358.75]]  

Multi-dimensional array of y's (lat)

[[-90.   -90.   -90.   ... -90.   -90.   -90.  ]  
...  
[ 90.    90.    90.   ...  90.    90.    90.  ]]

Multi-dimensional array of values

[[5076. 5076. 5076. ... 5076. 5076. 5076.]  
...  
[5138. 5138. 5138. ... 5138. 5138. 5138.]]

As these arrays are related one-to-one to each other, number of elements in each is 42048.

How do I create geoJSON from them?

In order to draw it chloropleth-like in Leaflet (here is more info on chloropleth-Leaflet if needed)

If I understand correctly, array of [x, y, val] will represent a Point and to draw the chloropleth (showing in colors values on the map) map I guess I need Polygons.

3
  • 1
    Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format. – PolyGeo Feb 16 '18 at 8:31
  • You are trying to convert a raster into a vector format, while it might be technically possible it is unlikely to be a good solution to your problem. – Ian Turton Feb 16 '18 at 9:15
  • @IanTurton, could you please advice in which direction to dig? – techkuz Feb 16 '18 at 9:37
1

Here is how to do it:

NOTE: You will need to install conda.

What this script does

  1. Takes in .zip full of grib2 and loops over them
  2. Makes some custom changes to the grib2 as a geopandas dataframe (you will need to customise this to your use case)
  3. Outputs in geojson with a .index file to make it easier for your front end to read it in.
import xarray as xr
import cfgrib
import eccodes
import geopandas as gpd
from shapely.geometry import Point
import json
import zipfile
import os
from os import walk

def extract_zip_file_and_return_list_of_file_names():
    cwd = os.getcwd()
    with zipfile.ZipFile(cwd + "/assets/YOURFILENAME.zip", 'r') as zip_ref:
        zip_ref.extractall(cwd + "/assets/")

    dir_path = cwd + "/assets/YOURDIRNAME/YOURFILENAME"

    list_of_file_names = []
    for (dirpath, dirnames, filenames) in walk(dir_path):
        list_of_file_names = filenames
        break

    return list_of_file_names



def process(list_of_grib2_filesnames):

    cwd = os.getcwd()
    list_of_json_filenames = []
    for filename in list_of_grib2_filesnames:
        print("Converting filename: " + filename)
        df = open_grib2_and_convert_to_dataframe(cwd + "/assets/YOURDIRNAME/" + filename)
        optomised_dataframe = optomise_dataframe(df)
        df_with_common_pressure_values = add_common_pressure_values_to_dataframe(optomised_dataframe)
        filename_in_json_format = write_to_json_file(df_with_common_pressure_values, filename)
        list_of_json_filenames.append(filename_in_json_format)

    write_to_json_file(list_of_json_filenames, cwd + "/index-file-name")

def open_grib2_and_convert_to_dataframe(file):
    data = xr.open_dataset(file, engine="cfgrib")
    # Convert grib2 into a geodataframe
    df = data.to_dataframe()
    return df


def optomise_dataframe(not_optomised_df):
    df = not_optomised_df.reset_index()
    try:
        geom = [Point(x,y) for x, y in zip(df['lon'], df['lat'])]
    except KeyError:
        geom = [Point(x,y) for x, y in zip(df['longitude'], df['latitude'])]

    gdf = gpd.GeoDataFrame(df, geometry=geom)


    del gdf['time']
    del gdf['step']
    del gdf['latitude']
    del gdf['longitude']
    del gdf['valid_time']

    # Delete empty values of turb
    optomised_dataframe = gdf[gdf.turb != 0.0]

    return optomised_dataframe


def _get_all_common_pressure_values_in_file(df):
    drop = df.drop_duplicates(subset=['isobaricInhPa'])
    common_pressure_values_in_df = []
    for pressure in drop['isobaricInhPa']:
        common_pressure_values_in_df.append(pressure)
    return common_pressure_values_in_df


def add_common_pressure_values_to_dataframe(df):
    geojson_data = {
        "geojson": json.loads(df.to_json()),
        "common": _get_all_common_pressure_values_in_file(df)
    }
    return geojson_data


def write_to_json_file(data, filename):
    with open(filename + ".json", 'w') as f:
        f.write(json.dumps(data))
        f.close()
    return filename + ".json"


if __name__ == "__main__":
    list_of_file_names = extract_zip_file_and_return_list_of_file_names()
    process(list_of_file_names)`

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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