4

I am trying to create a NetCDF from a CSV file. The idea comes from the fact that the original CSV comes with many columns, three of these representing:

  • a specific year ("time")
  • the Y coordinate of the observation in EPSG:4326 ("latitude")
  • the X coordinate of the observation in EPSG:4326 ("longitude")

As I want to be able to dynamically inspect each variable at a given time as a map, I want to convert it to NetCDF format.

I first import my CSV into a Pandas DataFrame, then I create a MultiIndex with the three dimensions I need, I convert this to an xarray DataSet, I write the CRS with rioxarray, and eventually I export it to a NetCDF nc file.

import pandas as pd
import rioxarray
import xarray as xr

df = pd.read_csv('my_data.csv')
df = df.set_index(['time', 'latitude', 'longitude'])
ds = df.to_xarray()
ds.rio.write_grid_mapping(inplace=True)
ds.rio.write_crs("epsg:4326", inplace=True)
ds.to_netcdf('my_data.nc')

However, when I try to import it as a mesh layer in QGIS, or in Panoply, the coordinates of my NetCDF are not seen correctly (see the pictures below).

I thought ds.rio.write_crs(4326, inplace=True) was enough to wtite the CRS to my dataset in order for it to be recognized by other software, but that does not seem enough.

I really wish to do it in Python to possibly avoid the use of other external software.

Info from ds.info:

Dimensions:           (latitude: 144, longitude: 140, time: 9)
Coordinates:
  * time              (time) datetime64[ns] 1980-01-01 1985-01-01 ... 2020-01-01
  * latitude          (latitude) float64 -34.75 -34.25 -33.75 ... 36.25 36.75
  * longitude         (longitude) float64 19.75 19.25 20.25 ... -23.75 -24.75
    spatial_ref       int64 0
Data variables:
    GHS_DENS_5avg     (time, latitude, longitude) float64 0.4566 nan ... nan nan
    rural             (time, latitude, longitude) float64 1.0 nan ... nan nan
    cities            (time, latitude, longitude) float64 0.0 nan ... nan nan
    temp_meanavg      (time, latitude, longitude) float64 17.4 nan ... nan nan
    temp_anom_5y      (time, latitude, longitude) float64 -0.01 nan ... nan nan
    prec_meanavg      (time, latitude, longitude) float64 37.1 nan ... nan nan
    prec_anom_5y      (time, latitude, longitude) float64 1.16 nan ... nan nan
    best              (time, latitude, longitude) float64 nan nan ... nan nan
    type_of_violence  (time, latitude, longitude) float64 nan nan ... nan nan
    n_event           (time, latitude, longitude) float64 nan nan ... nan nan
    any_event         (time, latitude, longitude) float64 0.0 nan ... nan nan>

Info from Panoply:

netcdf file:/home/umberto/Documents/jrc/teleworking/ca/projects/fabrizio/ciclim/data/Conflict%20data/africa_nc_conflict.nc {
  dimensions:
    time = 9;
    latitude = 144;
    longitude = 140;
  variables:
    long time(time=9);
      :units = "days since 1980-01-01 00:00:00";
      :calendar = "proleptic_gregorian";

    double latitude(latitude=144);
      :_FillValue = NaN; // double

    double longitude(longitude=140);
      :_FillValue = NaN; // double

    String iso3(time=9, latitude=144, longitude=140);
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double GHS_DENS_5avg(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :grid_mapping = "spatial_ref";
      :coordinates = "";

    double rural(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double cities(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double temp_meanavg(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :long_name = "5yr mean temperature";
      :units = "°C";
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double temp_anom_5y(time=9, latitude=144, longitude=140);
      :grid_mapping = "spatial_ref";
      :_FillValue = NaN; // double
      :coordinates = "";

    double prec_meanavg(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :grid_mapping = "spatial_ref";
      :coordinates = "";

    double prec_anom_5y(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double best(time=9, latitude=144, longitude=140);
      :grid_mapping = "spatial_ref";
      :coordinates = "";
      :_FillValue = NaN; // double
      :long_name = "Number of deaths";
      :units = "count";

    double type_of_violence(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :coordinates = "";
      :grid_mapping = "spatial_ref";

    double n_event(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :long_name = "Number of conflict events";
      :coordinates = "";
      :grid_mapping = "spatial_ref";
      :units = "count";

    double any_event(time=9, latitude=144, longitude=140);
      :_FillValue = NaN; // double
      :long_name = "Conflict events";
      :coordinates = "";
      :grid_mapping = "spatial_ref";
      :units = "Presence of conflict";

    long spatial_ref;
      :crs_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]";
      :semi_major_axis = 6378137.0; // double
      :semi_minor_axis = 6356752.314245179; // double
      :inverse_flattening = 298.257223563; // double
      :reference_ellipsoid_name = "WGS 84";
      :longitude_of_prime_meridian = 0.0; // double
      :prime_meridian_name = "Greenwich";
      :geographic_crs_name = "WGS 84";
      :grid_mapping_name = "latitude_longitude";
      :spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]";netcdf file:/home/umberto/Documents/jrc/teleworking/ca/projects/fabrizio/ciclim/data/Conflict%20data/africa_nc_conflict.nc {
  dimensions:
    time = 9;
    latitude = 144;
    longitude = 140;

  // global attributes:
}

NetCDF as mesh in QGIS

NetCDF as mesh in Panoply

1 Answer 1

1

Thanks to the answer of @snowman2, I was able to understand that the problem was that my coordinates were not properly sorted. Well, at least it seems to have solved the issue for Panoply and the plots I do with Python, while QGIS is still not displaying the layer in the correct position (but I guess this would be another question).

The updated code (with comment to emphasize the solution), is this one:

import pandas as pd
import rioxarray
import xarray as xr

df = pd.read_csv('my_data.csv')
df = df.set_index(['time', 'latitude', 'longitude'])
ds = df.to_xarray()
ds = ds.sortby(["time", "latitude", "longitude"]) # THIS SOLVED THE PROBLEM
ds.rio.write_crs("epsg:4326", inplace=True)
ds.to_netcdf('my_data.nc')

Panoply: enter image description here

QGIS (still bad, don't know why): enter image description here

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.