Getting rio.write_crs to work

[umbe1987]

Hello and first of all thanks for this neat project!

I am currently trying to create a NetCDF with Python and rioxarray that is in EPSG:4326 coordinate system.

This is a sort of alternative to my question already posted on gis.stackexchange.com.

I am starting from a CSV file with latitude, longitude, time columns which I use as dimensions, and other columns which are used as variables.

The workflow I am following to create the NetCDF file from the CSV is as follows:

• reading the CSV in a Pandas DataFrame
• setting the index of my df to be a MultiIndex using the three dimension columns (time, latitude, longitude, in this order)
• converting the df to xarray
• writing the grid_mapping with rio.write_grid_mapping
• writing the CRS with rio.write_crs
• exporting to NetCDF

Here is the code I am using:

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')

This all seems to write the nc file correctly, except that when I try to visualize it in QGIS or Panoply its content is not displayed in the correct place in the World.

It should cover Africa, while I am getting strange positioning like the ones in the pictures below. I would greatly appreciate any help to solve this issue. Many thanks in advance.

I also insert the info based on 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\"]]";

  // global attributes:
}

[snowman2]

Without the actual file, it is hard to really debug this as there are many potential issues.

1. Make sure they x/y coordinates are equally spaced and sorted after converting to xarray.
2. Make sure longitude values are from [-180, 180] and latitude are [-90, 90].

writing the grid_mapping with rio.write_grid_mapping

write_crs handles writing the grid mapping for you, so no need to do that in a separate step.

[umbe1987]

Thank you very much @snowman2 for your useful comment!

The data is equally spaced (0.5 degree), although I start from a grid like dataset which has some "holes" that I fill with this:

GRID_RESOLUTION = 0.5 # degree (CRS=EPSG:4326)

# create evenly-spaced coordinate pairs grid on the basis of lat lon values
x = np.arange(df['longitude'].min(), df['longitude'].max(), GRID_RESOLUTION)
y = np.arange(df['latitude'].min(), df['latitude'].max(), GRID_RESOLUTION)

df = df.set_index(['time', 'latitude', 'longitude'])

df = df.reindex(y, level='latitude')
df = df.reindex(x, level='longitude')

I think the problem might be in the sorting, as you suggested.

I normally sort my DataFrame with df.sort_index(inplace=True) before converting to xarray DataSet, but when I plot using python, I still have to do ds.sortby('longitude') otherwise I get the error:

The input coordinate is not sorted in increasing order along axis 0.

I also tried sorting directly the xarray DataSet using ds.sortby(["time", "latitude", "longitude"]), but I get the same error.

Do you know how am I supposed to sort my coordinates?

[umbe1987]

ok, sorry, apparently, I missed that ds.sortby() is not happening inplace, and as there's not an inplace argument, I had to reassing with ds = ds.sortby(["time", "latitude", "longitude"]), and the error is gone!

Thank you!

Still trying to figure out how to place the coordinates to let QGIS placing my data correctly but I guess I am close to the solution now, but at least PanoPly seems to like them ;)

[leomiquelutti]

hi @umbe1987,

did you solve it for QGIS? I am struggling with that for some days already.

My netcdf file does not show up in the correct location, while when I export my xarray as geotif with the to_raster method, it works!

[umbe1987]

hi @umbe1987,

did you solve it for QGIS? I am struggling with that for some days already.

My netcdf file does not show up in the correct location, while when I export my xarray as geotif with the to_raster method, it works!

@leomiquelutti I did not, but at least I am able to use the files within python just nicely. I used QGIS as a means to test whether what I was producing made sense, but apparently QGIS seemed to have an issue in representing (what most likely are) correct NetCDF files.