Converting NetCDF dataset with 2D coordinates to GeoTiff #47

goriliukasbuxton · 2019-09-17T17:59:26Z

import rioxarray
import xarray as xr

#Sentinel-5P data
xds = xr.open_dataset(r'S5P_NRTI_L2__NO2____20190513T181819_20190513T182319_08191_01_010301_20190513T185033.nc', group="/PRODUCT")

data = xds.nitrogendioxide_tropospheric_column

data.rio.to_raster("test1.tif")

gives error: "x dimension not found. 'set_spatial_dims()' can address this." rioxarray.exceptions.DimensionError: x dimension not found. 'set_spatial_dims()' can address this.

djhoese · 2019-09-17T18:05:08Z

What are the dimensions of data? As the error message says, it needs an 'x' dimension to know what to do. Your DataArray likely doesn't name its dimensions that way so you'll need to tell rioxarray which dimensions correspond to x and y by calling set_spatial_dims (data.rio.set_spatial_dims(...)).

goriliukasbuxton · 2019-09-17T18:13:39Z

this is the print out:
xarray.core.dataset.Dataset

<xarray.Dataset>
Dimensions:                                               (corner: 4, ground_pixel: 450, intensity_offset_polynomial_exponents: 1, layer: 34, polynomial_exponents: 6, scanline: 289, time: 1, vertices: 2)
Coordinates:
  * scanline                                              (scanline) float64 0.0 ... 288.0
  * ground_pixel                                          (ground_pixel) float64 0.0 ... 449.0
  * time                                                  (time) datetime64[ns] 2019-05-13
  * corner                                                (corner) float64 0.0 ... 3.0
  * polynomial_exponents                                  (polynomial_exponents) float64 0.0 ... 5.0
  * intensity_offset_polynomial_exponents                 (intensity_offset_polynomial_exponents) float64 0.0
  * layer                                                 (layer) float64 0.0 ... 33.0
  * vertices                                              (vertices) float64 0.0 1.0
    latitude                                              (time, scanline, ground_pixel) float32 ...
    longitude                                             (time, scanline, ground_pixel) float32 ...
Data variables:
    delta_time                                            (time, scanline) datetime64[ns] ...
    time_utc                                              (time, scanline) object ...
    qa_value                                              (time, scanline, ground_pixel) float32 ...
    nitrogendioxide_tropospheric_column                   (time, scanline, ground_pixel) float32 ...
    nitrogendioxide_tropospheric_column_precision         (time, scanline, ground_pixel) float32 ...
    nitrogendioxide_tropospheric_column_precision_kernel  (time, scanline, ground_pixel) float32 ...
    averaging_kernel                                      (time, scanline, ground_pixel, layer) float32 ...
    air_mass_factor_troposphere                           (time, scanline, ground_pixel) float32 ...
    air_mass_factor_total                                 (time, scanline, ground_pixel) float32 ...
    tm5_tropopause_layer_index                            (time, scanline, ground_pixel) float64 ...
    tm5_constant_a                                        (layer, vertices) float32 ...
    tm5_constant_b                                        (layer, vertices) float32 ...

djhoese · 2019-09-17T18:21:03Z

Here's the documentation for set_spatial_dims: https://corteva.github.io/rioxarray/html/rioxarray.html#rioxarray.rioxarray.XRasterBase.set_spatial_dims

snowman2 · 2019-09-17T18:36:28Z

Currently rioxarray assumes your spatial dimensions and coordinates are named the same.

For example:

<xarray.DataArray (band: 1, y: 245, x: 574)>
dask.array<shape=(1, 245, 574), dtype=float64, chunksize=(1, 245, 574)>
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 4.616e+06 4.616e+06 ... 4.615e+06 4.615e+06
  * x            (x) float64 4.25e+05 4.251e+05 ... 4.268e+05 4.268e+05
    spatial_ref  int64 0

@goriliukasbuxton, I see:

latitude (time, scanline, ground_pixel) float32 ...
longitude (time, scanline, ground_pixel) float32 ...

Can you define what scanline and ground_pixel are?

goriliukasbuxton · 2019-09-17T18:46:45Z

<xarray.DataArray 'ground_pixel' (ground_pixel: 450)>
array([  0.,   1.,   2., ..., 447., 448., 449.])
Coordinates:
  * ground_pixel  (ground_pixel) float64 0.0 1.0 2.0 3.0 ... 447.0 448.0 449.0
Attributes:
    units:      1
    axis:       X
    long_name:  across-track dimension index
    comment:    This coordinate variable defines the indices across track, fr...
<xarray.DataArray 'scanline' (scanline: 289)>
array([  0.,   1.,   2., ..., 286., 287., 288.])
Coordinates:
  * scanline  (scanline) float64 0.0 1.0 2.0 3.0 4.0 ... 285.0 286.0 287.0 288.0
Attributes:
    units:      1
    axis:       Y
    long_name:  along-track dimension index
    comment:    This coordinate variable defines the indices along track; ind...

snowman2 · 2019-09-17T19:09:45Z

Can you try: http://xarray.pydata.org/en/stable/generated/xarray.Dataset.rename_dims.html#xarray.Dataset.rename_dims

data = data.rename_dims({"ground_pixel": "longitude", "scanline": "latitude"})
data.rio.to_raster("test1.tif")

goriliukasbuxton · 2019-09-17T19:15:43Z

AttributeError: 'DataArray' object has no attribute 'rename_dims'

snowman2 · 2019-09-17T19:25:12Z

Whoops. wrong variable. How about this one:

xds = xds.rename_dims({"ground_pixel": "longitude", "scanline": "latitude"})
data = xds.nitrogendioxide_tropospheric_column
data.rio.to_raster("test1.tif")

goriliukasbuxton · 2019-09-17T19:29:58Z

no go:

site-packages\xarray\core\variable.py", line 1849, in __init__
    type(self).__name__)
ValueError: IndexVariable objects must be 1-dimensional

snowman2 · 2019-09-17T19:34:20Z

Probably naming conflicts. I wonder if you need to drop the ground_pixel and scanline coordinates first. Are you able to share the file?

goriliukasbuxton · 2019-09-17T19:38:11Z

https://drive.google.com/file/d/1tUStTSUoKoYpUAfzhDLwOSmebkYpp6Js/view?usp=sharing

snowman2 · 2019-09-17T19:58:24Z

Your data has unevenly spaced 2D lat/lon values.

That means that either:

Your data is not geographic and was re-projected to lat/lon in the 2D space to preserve the coordinate locations.
Your data is not represented in an evenly spaced grid.

If (1) is your problem, you will need to re-project your data from 2D latlon back to the projected coordinates in 1D, add the 1D x&y coordinates to your dataset. Then it should work.

If (2) is your problem, then you will need to decide on the resolution of your grid, the bounds of your area of interest, and resample your data to the new regularly spaced grid. Then it will work.

snowman2 · 2019-09-17T20:00:42Z

For (1), this may be helpful: https://github.com/snowman2/gazar/blob/7e0317466afa96ccdd5dd214b924b9c93d8a7788/gazar/grid.py#L683-L717

goriliukasbuxton · 2019-09-17T20:03:33Z

data from the same source, just would like to save it to files.
https://notes.stefanomattia.net/2018/03/27/Mount-Sinabung-eruption-as-seen-by-Sentinel-5P/
https://notes.stefanomattia.net/2019/05/21/Sentinel-5P-S-NPP-and-CH4-products-joint-plot/

leomiquelutti · 2020-08-05T12:13:17Z

Hi everyone, especially @djhoese, @snowman2 and @goriliukasbuxton

I am having similar problems. I use verde, from fatiando, for gridding and interpolating geospatial data.

Here is my example code. Example data can be downloaded here.

import pyproj
import verde as vd
import pandas
import rioxarray

# inputs
depth_of_interest = -273.4
spacing = 5e3
filename = 'F:\Google Drive Uff\Gawler Challenge\GC - 34\dados numericos\AusLAMP_MT_Gawler.xyzr'

# read data
data = pandas.read_csv(filename, delimiter=',')
data = data[data["depth_m"] == depth_of_interest]
print(data)

# re-project to UTM 
projection = pyproj.Proj(proj='utm',zone=53,ellps='WGS84', preserve_units=False)
proj_coords = projection(data.longitude.values, data.latitude.values)

# interpolate over regular grid
spline = vd.Spline().fit(proj_coords, data.log10_resistivity_ohm_m)
grid = spline.grid(spacing=spacing, data_names=["resistivity"])
grid = vd.distance_mask(proj_coords, maxdist=10e3, grid=grid)
print(grid)

# save to geotiff
grid.rio.set_crs("epsg:32753")
grid.rio.set_spatial_dims("easting", "northing", inplace=True)
# grid.rename({'easting': 'x','northing': 'y'})
grid["resistivity"].rio.to_raster('MT_25_verde.tif')

The error is always rioxarray.exceptions.DimensionError: x dimension not found. 'set_spatial_dims()' can address this., independtly whether I use:

grid.rio.set_spatial_dims("easting", "northing", inplace=True)
or grid.rename({'easting': 'x','northing': 'y'})

grid is an xarray, and the print(grid) outputs

<xarray.Dataset>
Dimensions:      (easting: 262, northing: 189)
Coordinates:
  * easting      (easting) float64 -1.279e+05 -1.229e+05 ... 1.172e+06 1.177e+06
  * northing     (northing) float64 -3.984e+06 -3.979e+06 ... -3.043e+06
Data variables:
    resistivity  (northing, easting) float64 -0.5191 -0.5184 nan ... 1.992 1.993
Attributes:
    metadata:  Generated by Spline(damping=None, engine='auto', force_coords=...

To check whether my data is over an evenly-spaced grid, I did the following checking and got these results:

grid.easting.values[1] - grid.easting.values[0] = 4998.4100816587015
grid.easting.values[100] - grid.easting.values[99] = 4998.4100816587015
grid.easting.values[180] - grid.easting.values[179] = 4998.410081658745

grid.northing.values[1] - grid.northing.values[0] = 5005.917018732522
grid.northing.values[100] - grid.northing.values[99] = 5005.917018732522
grid.northing.values[180] - grid.northing.values[179] = 5005.917018732056

My grid seems evenly spaced but for the very last decimals places. Any ideas on how to solve that, please?

snowman2 · 2020-08-05T13:02:48Z

@leomiquelutti what version of rioxarray are you using?

snowman2 · 2020-08-05T13:11:27Z

Also, it would be interesting to see the difference in interpolation between verde and geocube:

leomiquelutti · 2020-08-05T13:21:06Z

Hi @snowman2,

rioxarray version is 0.0.21. Should I update it?

Do you believe that the minor differences in the last decimals are responsible for that?

snowman2 · 2020-08-05T13:21:56Z

Should I update it?

Yes, please update to the latest and try it again.

snowman2 · 2020-08-05T13:22:34Z

Do you believe that the minor differences in the last decimals are responsible for that?

I don't believe so in this scenario.

snowman2 · 2020-08-05T13:32:16Z

import geopandas
import pandas
from geocube.api.core import make_geocube
from geocube.rasterize import rasterize_points_griddata

depth_of_interest = -273.4
spacing = 5e3
filename = 'AusLAMP_MT_Gawler.xyzr'

data = pandas.read_csv(filename, delimiter=',')
data = data[data["depth_m"] == depth_of_interest]

gdf = geopandas.GeoDataFrame(
    data, 
    geometry=geopandas.points_from_xy(data.longitude, data.latitude),
    crs="EPSG:4326"
)
gdf_utm = gdf.to_crs("EPSG:32753")

geo_grid = make_geocube(
    vector_data=gdf_utm,
    measurements=['log10_resistivity_ohm_m'],
    resolution=(-spacing, spacing),
    rasterize_function=rasterize_points_griddata,
)
geo_grid.log10_resistivity_ohm_m.rio.to_raster('MT_25_verde.tif')

snowman2 · 2020-08-05T13:34:33Z

$ geocube --show-versions
geocube v0.0.13

GDAL deps:
         fiona: 1.8.13
   GDAL[fiona]: 3.0.4
      rasterio: 1.1.5
GDAL[rasterio]: 3.0.4

Python deps:
       appdirs: 1.4.3
         click: 7.1.2
      datacube: 1.8.1
     geopandas: 0.8.0
     rioxarray: 0.0.31
        pyproj: 2.6.1.post1
        xarray: 0.15.1

System:
        python: 3.6.10 | packaged by conda-forge | (default, Apr 24 2020, 16:44:11)  [GCC 7.3.0]
    executable:~/miniconda/envs/geocube/bin/python3.6
       machine: Linux-4.15.0-112-generic-x86_64-with-debian-buster-sid

leomiquelutti · 2020-08-05T14:30:06Z

@snowman2 thanks for the support.

Unfortunately, updating rioxarray didn't solve my problem.
Your approach is very elegant. However, for several reasons, I need to do this with verde.

Any more ideas?

snowman2 · 2020-08-05T15:32:16Z

The issue you are running into is that the set_spatial_dims updated the dataset, but when you access the dataarray it no longer stores the variables from the dataset.

This works:

grid.resistivity.rio.set_spatial_dims("easting", "northing", inplace=True) \
    .rio.write_crs("EPSG:32753", inplace=True) \
    .rio.to_raster("MT_25_verde.tif")

Side note, you should probably do this to initialize Proj as it adds +south:

>>> import pyproj
>>> projection = pyproj.Proj("EPSG:32753")
>>> projection.srs
'+proj=utm +zone=53 +south +datum=WGS84 +units=m +no_defs'

snowman2 · 2020-08-05T15:38:58Z

Also, rename creates a copy and does not modify the original, so you need to do this for it to work:

grid = grid.rename({'easting': 'x','northing': 'y'})

leomiquelutti · 2020-08-06T21:52:46Z

This worked! The image appeared in the right location in ArcGIS! Thanks a lot, @snowman2!

I am still a beginner in Python, and would gladly appreciate some directions on the syntax you used (the slashes). A link would be very helpful, but don't bother if this requires lots of effort. I would like to understand this, to be clear.

grid.resistivity.rio.set_spatial_dims("easting", "northing", inplace=True) \
    .rio.write_crs("EPSG:32753", inplace=True) \
    .rio.to_raster("MT_25_verde.tif")

snowman2 · 2020-08-06T23:18:14Z

Glad to hear that it worked. Here is one reference: https://stackoverflow.com/a/53180/11622836

Bankimchandrayadav · 2020-12-07T14:00:02Z

@snowman2 Can you have a look at at this data: https://drive.google.com/file/d/1RoLPC6UaozsncO__qmst4sHrywN0_iyH/view?usp=sharing

I want to convert a single time band to geotiff in epsg:4326 projection.
For more ease please take a subset of 10x10 degrees like ds = ds.sel(lat=slice(25,35), lon=slice(75,85))

Please tell me how to convert it to tiff. I am having the same error about set_spatial_dims(). I tried setting x=lon and y=lat but not succeeding.

Bankimchandrayadav · 2020-12-08T16:24:57Z

@snowman2 Can you have a look at at this data: https://drive.google.com/file/d/1RoLPC6UaozsncO__qmst4sHrywN0_iyH/view?usp=sharing

I want to convert a single time band to geotiff in epsg:4326 projection.
For more ease please take a subset of 10x10 degrees like ds = ds.sel(lat=slice(25,35), lon=slice(75,85))

Please tell me how to convert it to tiff. I am having the same error about set_spatial_dims(). I tried setting x=lon and y=lat but not succeeding.

I believe I succeeded in converting it by changing the way I was opening the dataset:
Replacing:
ds = xarray.open_dataset('file.nc')
to
ds = rioxarray.open_rasterio('file.nc')
then typing these two lines:
ds = ds.rio.write_crs("epsg:4326") # crs set
ds['pre'].rio.to_raster('fileOutName.nc')

Please confirm if I did the right thing.
I have this trmm dataset where I am having a bit of confusion. I did convert it to .tiff by following this code:

ds = rioxarray.open_rasterio('trmmFile.nc')  
ds = ds.rename({'x': 'y','y': 'x'})  # renamed coz x and y were read interchanged
ds = ds.transpose('band', 'y', 'x')  # set dim order to as expected by rioxarray i.e. (band,y,x)
ds = ds.rio.write_crs("epsg:4326")  # crs set 
ds['precipitation'].rio.to_raster(outName) # saved as tif

TRMM file is here: https://drive.google.com/file/d/1g_Nf-KKNP1psIj3q9HdQHH6396yUBDsm/view?usp=sharing

snowman2 · 2020-12-08T16:49:07Z

This worked for me:

xds = xarray.open_dataset("1998_01_01.nc")
xds.precipitation.squeeze().transpose(
    'nlat', 'nlon'
).rio.set_spatial_dims(x_dim="nlon", y_dim="nlat").rio.to_raster("1998_01_01.tif")

Bankimchandrayadav · 2020-12-09T03:24:09Z

@snowman2 This worked. Thanks.

snowman2 · 2021-01-19T17:32:09Z

Closing for #209

danhamill mentioned this issue Mar 17, 2020

Export raster after reproject, not working. #96

Closed

snowman2 added the question Further information is requested label Apr 10, 2020

snowman2 changed the title ~~Converting NetCDF dataset array to GeoTiff using rioxarray, xarray Python~~ Converting NetCDF dataset with 2D coordinates to GeoTiff Oct 9, 2020

corteva locked as off-topic and limited conversation to collaborators Dec 22, 2020

snowman2 closed this as completed Jan 19, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Converting NetCDF dataset with 2D coordinates to GeoTiff #47

Converting NetCDF dataset with 2D coordinates to GeoTiff #47

goriliukasbuxton commented Sep 17, 2019 •

edited

Loading

djhoese commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019 •

edited by snowman2

Loading

djhoese commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020 •

edited

Loading

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020 •

edited

Loading

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 6, 2020

snowman2 commented Aug 6, 2020

Bankimchandrayadav commented Dec 7, 2020

Bankimchandrayadav commented Dec 8, 2020

snowman2 commented Dec 8, 2020

Bankimchandrayadav commented Dec 9, 2020

snowman2 commented Jan 19, 2021

Converting NetCDF dataset with 2D coordinates to GeoTiff #47

Converting NetCDF dataset with 2D coordinates to GeoTiff #47

Comments

goriliukasbuxton commented Sep 17, 2019 • edited Loading

djhoese commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019 • edited by snowman2 Loading

djhoese commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

snowman2 commented Sep 17, 2019

snowman2 commented Sep 17, 2019

goriliukasbuxton commented Sep 17, 2019

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020

snowman2 commented Aug 5, 2020 • edited Loading

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 5, 2020

snowman2 commented Aug 5, 2020 • edited Loading

snowman2 commented Aug 5, 2020

leomiquelutti commented Aug 6, 2020

snowman2 commented Aug 6, 2020

Bankimchandrayadav commented Dec 7, 2020

Bankimchandrayadav commented Dec 8, 2020

snowman2 commented Dec 8, 2020

Bankimchandrayadav commented Dec 9, 2020

snowman2 commented Jan 19, 2021

goriliukasbuxton commented Sep 17, 2019 •

edited

Loading

goriliukasbuxton commented Sep 17, 2019 •

edited by snowman2

Loading

snowman2 commented Aug 5, 2020 •

edited

Loading

snowman2 commented Aug 5, 2020 •

edited

Loading