Clipped netcdf produces nan values #672

gkb999 · 2023-05-07T08:29:45Z

gkb999
May 7, 2023

Greetings,
Im using rioxarray to clip my netcdf file based on shapefile.
Although I had issues initially, I could successfully clip it.
However, when I'm exporting it as a netcdf file and then check the netcdf file in other sofwtare, it shows all nan values.
When I try opening within python, I can plot, and below is an example of the mask.

The values within (in panoply)

What could be the issue?
Mask Code

poly = regionmask.Regions(np.array(sf1.geometry)) # my shape
mask = poly.mask(ds['z'].isel(time=0),lat_name = 'y',lon_name = 'x', wrap_lon=False) #masked

While clipping, I did not use any like for masking nan values. Is there anything as such?

When I look into "mask":

mask

<xarray.DataArray 'mask' (y: 1328, x: 1264)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
  * y        (y) float64 -3.947e+06 -3.941e+06 ... 4.341e+06 4.347e+06
  * x        (x) float64 -3.947e+06 -3.941e+06 ... 3.941e+06 3.947e+06

How to get rid of those empty spaces
How to convert x and y to lat-lon format.

I tried using transform utility form pyproj - that was a mess and even masking wasn't successful.

Thanks in advance
Warm Regards

Answered by gkb999

May 18, 2023

Thanks for getting back @ccalvocm
The original data's

when i type :
ds.polar_stereographic.attrs

{'grid_mapping_name': 'polar_stereographic',
 'straight_vertical_longitude_from_pole': 0.0,
 'false_easting': 0.0,
 'false_northing': 0.0,
 'latitude_of_projection_origin': -90.0,
 'standard_parallel': -70.0,
 'long_name': 'CRS definition',
 'longitude_of_prime_meridian': 0.0,
 'semi_major_axis': 6378273.0,
 'inverse_flattening': 298.279411123064,
 'spatial_ref': 'PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7…

View full answer

ccalvocm · 2023-05-17T00:49:18Z

ccalvocm
May 17, 2023

Greetings @gkb999,

Maybe the data shows a lot of NaNs, but not all of the array data is NaN. ¿Can you check that with?:

import numpy as np
np.nanmax(YourXArray.values)

Maybe the software you are using to open the exported NetCDF overlooks few not NaN values.

Another procedure for clipping 100s of nc files is CDO:

Here's CDO:
https://code.mpimet.mpg.de/projects/cdo/embedded/index.html

Here's how to mask NetCDFs with CDO:
Mask NETCDF file by shapefile?
Masking a netcdf file by CDO

Video:
Video

Kind regards

0 replies

snowman2 · 2023-05-17T01:45:00Z

snowman2
May 17, 2023
Maintainer

@gkb999, are you able to provide a simple re-producible code snippet with input files?

ANSWER:

import xarray
import rioxarray
import geopandas

xds = xarray.open_mfdataset('H:/*.nc',combine = 'nested', concat_dim="time")
xds = xds.rio.write_crs("EPSG:3412", inplace=True)
geodf = geopandas.read_file('H:/roi_clip/ROI.shp')
clip = xds.rio.clip(geodf.geometry.values, geodf.crs).rio.reproject("EPSG:4326")

11 replies

ccalvocm May 18, 2023

Greetings @gkb999

To convert to lat-lon you would need to reproject your xarray:

import numpy as np
from pyproj import Transformer
xv, yv = np.meshgrid(sic_daily.x, sic_daily.y)
transformer = Transformer.from_crs("epsg:YOUREPSGNUMBERHERE", "epsg:4326", always_xy=True, )
lon, lat = transformer.transform(xv, yv)
sic_daily.coords['lon'] = (('y', 'x'), lon)
sic_daily.coords['lat'] = (('y', 'x'), lat)
sic_daily.attrs['crs'] = '+init=epsg:4326'

gkb999 May 18, 2023
Author

Many thanks for the prompt reply @ccalvocm
It definetely convereted the x-y to lat-lon, but the converesion appears wwrong to me.
What I did:

from pyproj import Transformer
xv, yv = np.meshgrid(sic_daily.x, sic_daily.y)
transformer = Transformer.from_crs("epsg:3412",  "epsg:4326", always_xy=True, ) 
lon, lat = transformer.transform(xv, yv)
sic_daily.coords['lon'] = (('y', 'x'), lon)
sic_daily.coords['lat'] = (('y', 'x'), lat)
sic_daily.attrs['crs']  = '+init=epsg:4326'

My ROI's spatial Reference which is same as the raster clipped:

The Output:

<xarray.Dataset>
Dimensions:    (x: 1264, y: 1328, time: 3287, bnds: 2)
Coordinates:
  * x          (x) float64 -3.947e+06 -3.941e+06 ... 3.941e+06 3.947e+06
  * y          (y) float64 -3.947e+06 -3.941e+06 ... 4.341e+06 4.347e+06
  * time       (time) int32 0 1 2 3 4 5 6 ... 3280 3281 3282 3283 3284 3285 3286
    lon        (y, x) float64 -135.0 -135.0 -135.1 -135.1 ... 42.15 42.19 42.24
    lat        (y, x) float64 -41.48 -41.51 -41.55 ... -39.33 -39.3 -39.26
Dimensions without coordinates: bnds
Data variables:
    z          (time, y, x) float32 ...
    x_bnds     (x, bnds) float64 -3.95e+06 -3.944e+06 ... 3.944e+06 3.95e+06
    time_bnds  (time, bnds) float64 -0.5 0.5 0.5 ... 3.286e+03 3.286e+03
Attributes:
    crs:      +init=epsg:4326

The lat-lon should not exceed 200
the shape bounds are typically : longitude: 160E to 200 and latitude : 90S to -75S
Also, the map when I plot:

I need something with lat-lon values like:

ccalvocm May 18, 2023

Greetings @gkb999 ,

¿Are you very sure your original data is EPSG:3412?, I can see values of x ranging to 3.947e+06 and y to 4.347e+06

Kind regards

gkb999 May 18, 2023
Author

Thanks for getting back @ccalvocm
The original data's

when i type :
ds.polar_stereographic.attrs

{'grid_mapping_name': 'polar_stereographic',
 'straight_vertical_longitude_from_pole': 0.0,
 'false_easting': 0.0,
 'false_northing': 0.0,
 'latitude_of_projection_origin': -90.0,
 'standard_parallel': -70.0,
 'long_name': 'CRS definition',
 'longitude_of_prime_meridian': 0.0,
 'semi_major_axis': 6378273.0,
 'inverse_flattening': 298.279411123064,
 'spatial_ref': 'PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","6054"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4054"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3412"]]',
 'GeoTransform': '-3950000 6250 0 4350000 0 -6250 '}

Answer selected by gkb999

gkb999 May 18, 2023
Author

Can probably be anti-meridian issues?
The data crosses 180E, so could that be?
But, either way, even latitude-values aren't in 40S-50S range (originally : 90S-72S)

snowman2 May 18, 2023
Maintainer

I don't see any code using rioxarray yet. Did you follow these examples?

gkb999 May 19, 2023
Author

Hi, the clip_geom shows the example of raster input (tiff format).
I'll be a bit clear.

My data is in netcdf format

<xarray.Dataset>
Dimensions:              (time: 30, x: 1264, y: 1328)
Coordinates:
  * x                    (x) float64 -3.947e+06 -3.941e+06 ... 3.947e+06
  * y                    (y) float64 -3.947e+06 -3.941e+06 ... 4.347e+06
Dimensions without coordinates: time
Data variables:
    polar_stereographic  (time) |S1 b'' b'' b'' b'' b'' ... b'' b'' b'' b'' b''
    z                    (time, y, x) float32 dask.array<chunksize=(1, 1328, 1264), meta=np.ndarray>
Attributes:
    author:       University of Bremen, Gunnar Spreen [gunnar.spreen@uni-brem...
    Conventions:  CF-1.5
    GMT_version:  5.2.1 (r15220) [64-bit] [MP]
    history:      Mon Jan 21 20:08:59 2019: GDAL CreateCopy( /ssmi/www/htdocs...
    info:         AMSR2 sea ice concentration based on the ASI algorithm (Spr...
    title:        Produced by grdmath
    GDAL:         GDAL 2.1.3, released 2017/20/01

Original projection info:

{'grid_mapping_name': 'polar_stereographic',
 'straight_vertical_longitude_from_pole': 0.0,
 'false_easting': 0.0,
 'false_northing': 0.0,
 'latitude_of_projection_origin': -90.0,
 'standard_parallel': -70.0,
 'long_name': 'CRS definition',
 'longitude_of_prime_meridian': 0.0,
 'semi_major_axis': 6378273.0,
 'inverse_flattening': 298.279411123064,
 'spatial_ref': 'PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","6054"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4054"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3412"]]',
 'GeoTransform': '-3950000 6250 0 4350000 0 -6250 '}

For clipping this data with my shapefile (ROI)

import xarray
import rioxarray
import geopandas
from shapely.geometry import mapping
###############
xds = xr.open_mfdataset('H:/*.nc',combine = 'nested', concat_dim="time")
xds = xds.rio.write_crs("EPSG:3412", inplace=True)
geodf = geopandas.read_file('H:/roi_clip/ROI.shp')
####
geometries = [
    {
        'type': 'Polygon',
        'coordinates': [[-510168.6036424115300179,1119413.0793986767530441], [515779.1901665925979614,1888217.9479051679372787]]
    }
]
geometries = geodf.geometry.apply(mapping)
clip = xds.rio.clip(geodf.geometry.values, geodf.crs, drop=False, invert=True)

which doesn't work

xarray.Dataset
Dimensions:
x: 1264y: 1328time: 30
Coordinates:
x
(x)
float64
-3.947e+06 -3.941e+06 ... 3.947e+06
y
(y)
float64
-3.947e+06 -3.941e+06 ... 4.347e+06
polar_stereographic
()
int32
0
Data variables:
z
(time, y, x)
float32
dask.array<chunksize=(1, 1328, 1264), meta=np.ndarray>
Indexes: (2)
Attributes:
author :
University of Bremen, Gunnar Spreen [gunnar.spreen@uni-bremen.de], Christian Melsheimer [melsheim@uni-bremen.de]
Conventions :
CF-1.5
GMT_version :
5.2.1 (r15220) [64-bit] [MP]
history :
Mon Jan 21 20:08:59 2019: GDAL CreateCopy( /ssmi/www/htdocs/data/amsr2/asi_daygrid_swath/s6250/netcdf/2016/asi-AMSR2-s6250-20160401-v5.4.nc, ... )
grdmath ../s6250/2016/apr/Antarctic/asi-AMSR2-s6250-20160401-v5.4.hdf FLIPUD = tmp1.nc
info :
AMSR2 sea ice concentration based on the ASI algorithm (Spreen et al., 2008).
title :
Produced by grdmath
GDAL :
GDAL 2.1.3, released 2017/20/01

and looks just the same:
#############
Now I set the crs and tried:

ds.rio.set_crs("epsg:3412")
geodf.geometry
clipped = ds.rio.clip(sf1.geometry.apply(mapping), sf1.crs, drop=False)

it worked

But, those blanks or those empty spaces

I need a zoomed in map (first priority)
Second, Convert those meters into lat-lon (Which constantly shows successful but is wrong.

i TRIED:

from pyproj import Transformer
xv, yv = np.meshgrid(sic_daily.x, sic_daily.y)
transformer = Transformer.from_crs("epsg:3412",  "epsg:4326", always_xy=True) 
lon, lat = transformer.transform(xv, yv)
sic_daily.coords['lon'] = (('y', 'x'), lon)
sic_daily.coords['lat'] = (('y', 'x'), lat)
sic_daily.attrs['crs']  = '+init=epsg:4326'

Lat-lon converted are wrong.

<xarray.Dataset>
Dimensions:    (x: 1264, y: 1328, time: 3287, bnds: 2)
Coordinates:
  * x          (x) float64 -3.947e+06 -3.941e+06 ... 3.941e+06 3.947e+06
  * y          (y) float64 -3.947e+06 -3.941e+06 ... 4.341e+06 4.347e+06
  * time       (time) int32 0 1 2 3 4 5 6 ... 3280 3281 3282 3283 3284 3285 3286
    lon        (y, x) float64 -135.0 -135.0 -135.1 -135.1 ... 42.15 42.19 42.24
    lat        (y, x) float64 -41.48 -41.51 -41.55 ... -39.33 -39.3 -39.26
Dimensions without coordinates: bnds
Data variables:
    z          (time, y, x) float32 ...
    x_bnds     (x, bnds) float64 -3.95e+06 -3.944e+06 ... 3.944e+06 3.95e+06
    time_bnds  (time, bnds) float64 -0.5 0.5 0.5 ... 3.286e+03 3.286e+03
Attributes:
    crs:      +init=epsg:4326

snowman2 May 19, 2023
Maintainer

If you set drop=False, it acts like a mask and doesn't clip to the geometry. What you want is:

clipped = ds.rio.clip(sf1.geometry.values, sf1.crs)

gkb999 · 2023-05-19T04:33:23Z

gkb999
May 19, 2023
Author

Just one last question:
I'm unable to export this to netcdf:

rs_clip_reproj.to_netcdf('H:/AMSR2_rs_latlon.nc')

gives the following error:

ValueError                                Traceback (most recent call last)
Cell In[92], line 2
      1 #rs_clip.to_netcdf('H:/AMSR2_rs_m.nc')
----> 2 rs_clip_reproj.to_netcdf('H:/AMSR2_rs_latlon.nc')

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\core\dataarray.py:3952, in DataArray.to_netcdf(self, path, mode, format, group, engine, encoding, unlimited_dims, compute, invalid_netcdf)
   3948 else:
   3949     # No problems with the name - so we're fine!
   3950     dataset = self.to_dataset()
-> 3952 return to_netcdf(  # type: ignore  # mypy cannot resolve the overloads:(
   3953     dataset,
   3954     path,
   3955     mode=mode,
   3956     format=format,
   3957     group=group,
   3958     engine=engine,
   3959     encoding=encoding,
   3960     unlimited_dims=unlimited_dims,
   3961     compute=compute,
   3962     multifile=False,
   3963     invalid_netcdf=invalid_netcdf,
   3964 )

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\backends\api.py:1217, in to_netcdf(dataset, path_or_file, mode, format, group, engine, encoding, unlimited_dims, compute, multifile, invalid_netcdf)
   1212 # TODO: figure out how to refactor this logic (here and in save_mfdataset)
   1213 # to avoid this mess of conditionals
   1214 try:
   1215     # TODO: allow this work (setting up the file for writing array data)
   1216     # to be parallelized with dask
-> 1217     dump_to_store(
   1218         dataset, store, writer, encoding=encoding, unlimited_dims=unlimited_dims
   1219     )
   1220     if autoclose:
   1221         store.close()

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\backends\api.py:1264, in dump_to_store(dataset, store, writer, encoder, encoding, unlimited_dims)
   1261 if encoder:
   1262     variables, attrs = encoder(variables, attrs)
-> 1264 store.store(variables, attrs, check_encoding, writer, unlimited_dims=unlimited_dims)

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\backends\common.py:267, in AbstractWritableDataStore.store(self, variables, attributes, check_encoding_set, writer, unlimited_dims)
    264 if writer is None:
    265     writer = ArrayWriter()
--> 267 variables, attributes = self.encode(variables, attributes)
    269 self.set_attributes(attributes)
    270 self.set_dimensions(variables, unlimited_dims=unlimited_dims)

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\backends\common.py:356, in WritableCFDataStore.encode(self, variables, attributes)
    353 def encode(self, variables, attributes):
    354     # All NetCDF files get CF encoded by default, without this attempting
    355     # to write times, for example, would fail.
--> 356     variables, attributes = cf_encoder(variables, attributes)
    357     variables = {k: self.encode_variable(v) for k, v in variables.items()}
    358     attributes = {k: self.encode_attribute(v) for k, v in attributes.items()}

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\conventions.py:891, in cf_encoder(variables, attributes)
    888 # add encoding for time bounds variables if present.
    889 _update_bounds_encoding(variables)
--> 891 new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
    893 # Remove attrs from bounds variables (issue #2921)
    894 for var in new_vars.values():

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\conventions.py:891, in <dictcomp>(.0)
    888 # add encoding for time bounds variables if present.
    889 _update_bounds_encoding(variables)
--> 891 new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
    893 # Remove attrs from bounds variables (issue #2921)
    894 for var in new_vars.values():

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\conventions.py:304, in encode_cf_variable(var, needs_copy, name)
    301 var = ensure_dtype_not_object(var, name=name)
    303 for attr_name in CF_RELATED_DATA:
--> 304     pop_to(var.encoding, var.attrs, attr_name)
    305 return var

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\coding\variables.py:134, in pop_to(source, dest, key, name)
    132 value = source.pop(key, None)
    133 if value is not None:
--> 134     safe_setitem(dest, key, value, name=name)
    135 return value

File ~\Anaconda3\envs\xMCA\lib\site-packages\xarray\coding\variables.py:115, in safe_setitem(dest, key, value, name)
    113 if key in dest:
    114     var_str = f" on variable {name!r}" if name else ""
--> 115     raise ValueError(
    116         "failed to prevent overwriting existing key {} in attrs{}. "
    117         "This is probably an encoding field used by xarray to describe "
    118         "how a variable is serialized. To proceed, remove this key from "
    119         "the variable's attributes manually.".format(key, var_str)
    120     )
    121 dest[key] = value

ValueError: failed to prevent overwriting existing key grid_mapping in attrs. This is probably an encoding field used by xarray to describe how a variable is serialized. To proceed, remove this key from the variable's attributes manually.

After I reproject, I wrapped the longitudes to 0-360 degrees, and it shows all Nan and the image is like:

Before 360 conversion (it has a few values)

<xarray.DataArray 'z' (time: 3287, y: 5, x: 240)>
array([[[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       , 29.743046 , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [80.83667  , 96.64541  ,  0.       , ...,        nan,
                nan, 89.27075  ]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       , 34.537308 , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [89.776985 , 94.760506 ,  1.0094886, ...,        nan,
...
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [90.618095 , 90.65166  ,  0.       , ...,        nan,
                nan, 92.54029  ]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [84.26544  , 90.546814 ,  0.8791098, ...,        nan,
                nan, 86.60795  ]]], dtype=float32)
Coordinates:
  * x                    (x) float64 -179.3 -177.8 -176.3 ... 175.9 177.4 178.9
  * y                    (y) float64 -72.12 -73.62 -75.12 -76.62 -78.12

After wrapping lon to 0-360 range, its all nan

<xarray.DataArray 'z' (time: 3287, y: 5, x: 240)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       ...,

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * x                    (x) float64 0.5872 2.086 3.585 ... 356.1 357.6 359.1
  * y                    (y) float64 -72.12 -73.62 -75.12 -76.62 -78.12

Also, I generally use the Orthographic projection for plotting these data.
But, although the environment has cartopy it still says

AttributeError                            Traceback (most recent call last)
Cell In[30], line 4
      2 from cartopy.feature import LAND
      3 #import cartopy as ccrs
----> 4 p = latlon.isel(time=21).plot(subplot_kws=dict(projection=ccrs.Orthographic()),transform=ccrs.PlateCarree())
      5 p.axes.coastlines()

AttributeError: module 'cartopy' has no attribute 'Orthographic'

3 replies

gkb999 May 19, 2023
Author

Well, I got why the image after reprojection was so terrible @snowman2
the longitudes were not reprojected right.
As the data passes by 180 (anti-meridian issues? as reprojection was in 4326?)
If you see this

<xarray.DataArray 'z' (time: 3287, y: 5, x: 240)>
array([[[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       , 29.743046 , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [80.83667  , 96.64541  ,  0.       , ...,        nan,
                nan, 89.27075  ]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       , 34.537308 , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [89.776985 , 94.760506 ,  1.0094886, ...,        nan,
...
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [90.618095 , 90.65166  ,  0.       , ...,        nan,
                nan, 92.54029  ]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [ 0.       ,        nan,        nan, ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [ 0.       ,  0.       ,  0.       , ...,  0.       ,
          0.       ,  0.       ],
        [84.26544  , 90.546814 ,  0.8791098, ...,        nan,
                nan, 86.60795  ]]], dtype=float32)
Coordinates:
  * x                    (x) float64 -179.3 -177.8 -176.3 ... 175.9 177.4 178.9
  * y                    (y) float64 -72.12 -73.62 -75.12 -76.62 -78.12

the x coordinate is showing -179 to +179. That means, the data there was never clipped.

snowman2 May 24, 2023
Maintainer

Related: #230

gkb999 May 24, 2023
Author

Hi @snowman2 thanks for the related content.
I'm anyway working these clips and reprojections.
I realized, the clip worked the other day but the longitude values are not right. When I used polar-projection, it appeared exactly where it is supposed to be.

but the longitudes show 0-360 (as I applied the wrap)

How to get the right longitude values there?
I thought reprojecting before clipping might sustain the longitude values.

So, what I did was, using the same method as mentioned by you,

ds.rio.set_crs("EPSG:3412", inplace=True) #set crs 
ds.z.rio.reproject("EPSG:4326") # reprojected in lat-lon

However this only works if I'm passing ds.z (that is data-array) and not on dataset directly

ds.rio.reproject("EPSG:4326")

###gives following error
MissingSpatialDimensionError              Traceback (most recent call last)
Cell In[258], line 1
----> 1 ds_latlon = ds_meters.rio.reproject("EPSG:4326")

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\rioxarray\raster_dataset.py:118, in RasterDataset.reproject(self, dst_crs, resolution, shape, transform, resampling, nodata, **kwargs)
    114 try:
    115     x_dim, y_dim = _get_spatial_dims(self._obj, var)
    116     resampled_dataset[var] = (
    117         self._obj[var]
--> 118         .rio.set_spatial_dims(x_dim=x_dim, y_dim=y_dim, inplace=True)
    119         .rio.reproject(
    120             dst_crs,
    121             resolution=resolution,
    122             shape=shape,
    123             transform=transform,
    124             resampling=resampling,
    125             nodata=nodata,
    126             **kwargs,
    127         )
    128     )
    129 except MissingSpatialDimensionError:
    130     if len(self._obj[var].dims) >= 2 and not get_option(
    131         SKIP_MISSING_SPATIAL_DIMS
    132     ):

File ~\Anaconda3\envs\xcdat_env\Lib\site-packages\rioxarray\rioxarray.py:864, in XRasterBase.set_spatial_dims(self, x_dim, y_dim, inplace)
    862     data_obj.rio._y_dim = y_dim
    863 else:
--> 864     raise MissingSpatialDimensionError(
    865         f"y dimension ({y_dim}) not found.{_get_data_var_message(data_obj)}"
    866     )
    867 return data_obj

MissingSpatialDimensionError: y dimension (y) not found. Data variable: x_bnds

I converetd the data-array to dataset

ds_latlon = ds_latlon.to_dataset()

but when I save, I don't see any array (it's just a single array) - no time - dimension left which I added manually before saving.
I'm using

del ds.z.attrs['grid_mapping'] #deleting the attribute as it is causing error while saving which is commonly reported across stack-flow and found deleting attribute is a solution
ds.to_netcdf('H:/latlon.nc')

Also, how do I remove nan-values from the final product?

My questions are:

Why am i unable to reproject dataset dircetly and only data-ararys?
How to completely not use nan-values for later analysis
Why only single dimension is being saved?

Many thanks in advance

snowman2 · 2023-05-24T13:21:34Z

snowman2
May 24, 2023
Maintainer

Why am i unable to reproject dataset dircetly and only data-ararys?

You can sometimes re-project Datasets. rioxarray tries to guess which DataArray's within a dataset are able to be re-projected. Sometimes, it thinks there is a bad array (in your scenario). If you only want the geospatial data, you can pre-select it from the dataset before re-projecting.

How to completely not use nan-values for later analysis

Mind re-wording/clarifying this question. Not sure what you mean.

Why only single dimension is being saved?

Not sure. The time dimension should be preserved ... hard to debug without the input files.

1 reply

gkb999 May 24, 2023
Author

You can sometimes re-project Datasets. rioxarray tries to guess which DataArray's within a dataset are able to be re-projected. Sometimes, it thinks there is a bad array (in your scenario). If you only want the geospatial data, you can pre-select it from the dataset before re-projecting.

Okay, I did the way mentioned, pre-selected/clipped before reprojection, but the longitudes are arranged as -180 to +180 , although the data doesn't contain those longitudes.
Can this be sorted? Or should be rearrange / manipulate longitude coordinates manually?

Mind re-wording/clarifying this question. Not sure what you mean.

Sorry for not making sense initially. I wanted to ask that when I'm saving back the data clipped and reprojected, a lot of nan values come along. So, is there a way while clipping to ignore these nans?

Anyway, that isn't of concern too - as I can use xarray to get rid of NaNs.

Not sure. The time dimension should be preserved ... hard to debug without the input files.

The time-dimensions is sorted. So, this part isn't an issue anymore

The only issue that's bothering is, sorting the longitudes coordinates.

gkb999 · 2023-05-25T01:50:23Z

gkb999
May 25, 2023
Author

Hi again,
Is it that the reprojection would distort and change the data resolution?

Is there a way we can just use the existing x-y grid is meters to generate lat-lon and have additional coordinates there?

Similar Links: https://gis.stackexchange.com/questions/245151/why-does-reprojection-with-gdalwarp-change-pixel-size
that talk about distortions at poles with wgs84.

5 replies

snowman2 May 25, 2023
Maintainer

You can pass in resolution= to reproject.

gkb999 May 25, 2023
Author

You can pass in resolution= to reproject.

rep = ds.rio.reproject(
    "EPSG:4326",
    resolution=6.25, #in kilometers
    nodata=np.nan # Set NaN as NoData
)

Is this right?
This doesn't seem to work either

Can I use

reproject_match 
``` ???

snowman2 May 25, 2023
Maintainer

The resolution is the resolution in the destination CRS. So, the units are degrees.

gkb999 May 25, 2023
Author

Thanks for the option @snowman2 but maybe this doesn't work
Firstly, this is never success because of units conversion - as the meters is a distance and degrees is angle, it changes as it goes towards the poles, especially the longitude. So, should accept that rioxarray may not be able to just by passing resolution
Secondly, it runs into memory issue, although the file isn't huge.
Gridded data is obviously higher in space.
Although I found a function that can convert xy to polar coordinates (https://github.com/nsidc/polarstereo-lonlat-convert-py/blob/main/README.md), but again, the coordinate are wrong there.

Anyway, thanks for the assistance so far, and if the resolution option works, please show me some sample (mine is 6.25km for my data).

snowman2 May 25, 2023
Maintainer

Firstly, this is never success because of units conversion - as the meters is a distance and degrees is angle, it changes as it goes towards the poles, especially the longitude. So, should accept that rioxarray may not be able to just by passing resolution

The grid will be equally spaced in the projection you choose. So, it will need equal intervals in degrees. As you have concluded, this will not be equal in meters. So, you will have to chose a value that is the closest approximation for what you want or pick an equal distance projection.

Clipped netcdf produces nan values #672

gkb999 May 7, 2023

Replies: 5 comments · 20 replies

ccalvocm May 17, 2023

snowman2 May 17, 2023 Maintainer

ccalvocm May 18, 2023

gkb999 May 18, 2023 Author

ccalvocm May 18, 2023

gkb999 May 18, 2023 Author

gkb999 May 18, 2023 Author

snowman2 May 18, 2023 Maintainer

gkb999 May 19, 2023 Author

snowman2 May 19, 2023 Maintainer

gkb999 May 19, 2023 Author

gkb999 May 19, 2023 Author

snowman2 May 24, 2023 Maintainer

gkb999 May 24, 2023 Author

snowman2 May 24, 2023 Maintainer

gkb999 May 24, 2023 Author

gkb999 May 25, 2023 Author

snowman2 May 25, 2023 Maintainer

gkb999 May 25, 2023 Author

snowman2 May 25, 2023 Maintainer

gkb999 May 25, 2023 Author

snowman2 May 25, 2023 Maintainer

gkb999
May 7, 2023

Replies: 5 comments 20 replies

ccalvocm
May 17, 2023

snowman2
May 17, 2023
Maintainer

gkb999 May 18, 2023
Author

gkb999 May 18, 2023
Author

gkb999 May 18, 2023
Author

snowman2 May 18, 2023
Maintainer

gkb999 May 19, 2023
Author

snowman2 May 19, 2023
Maintainer

gkb999
May 19, 2023
Author

gkb999 May 19, 2023
Author

snowman2 May 24, 2023
Maintainer

gkb999 May 24, 2023
Author

snowman2
May 24, 2023
Maintainer

gkb999 May 24, 2023
Author

gkb999
May 25, 2023
Author

snowman2 May 25, 2023
Maintainer

gkb999 May 25, 2023
Author

snowman2 May 25, 2023
Maintainer

gkb999 May 25, 2023
Author

snowman2 May 25, 2023
Maintainer