5

I have a NetCDF .nc file that gives me an initial error Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute from gdalinfo.

After searching I found out this relates to NetCDF files having no inherent CRS. Whats the best way to add a CRS, e.g. EPSG:4326, etc.? I've seen the rioxarray package mentioned but I had some trouble with the to_raster function.

Here is the data:

$ gdalinfo Flood.nc
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: Flood.nc
Size is 512, 512
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"Flood.nc":Pressure
  SUBDATASET_1_DESC=[155x155] Pressure (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"Flood.nc":Radar
  SUBDATASET_2_DESC=[155x155] Radar (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"Flood.nc":Rain6h
  SUBDATASET_3_DESC=[155x155] Rain6h (32-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"Flood.nc":RainTot
  SUBDATASET_4_DESC=[155x155] RainTot (32-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"Flood.nc":RH
  SUBDATASET_5_DESC=[155x155] RH (32-bit floating-point)
  SUBDATASET_6_NAME=NETCDF:"Flood.nc":U
  SUBDATASET_6_DESC=[155x155] U (32-bit floating-point)
  SUBDATASET_7_NAME=NETCDF:"Flood.nc":V
  SUBDATASET_7_DESC=[155x155] V (32-bit floating-point)
  SUBDATASET_8_NAME=NETCDF:"Flood.nc":Latitude
  SUBDATASET_8_DESC=[155x155] Latitude (32-bit floating-point)
  SUBDATASET_9_NAME=NETCDF:"Flood.nc":Longitude
  SUBDATASET_9_DESC=[155x155] Longitude (32-bit floating-point)
  SUBDATASET_10_NAME=NETCDF:"Flood.nc":Terrain
  SUBDATASET_10_DESC=[155x155] Terrain (32-bit floating-point)
  SUBDATASET_11_NAME=NETCDF:"Flood.nc":Landuse
  SUBDATASET_11_DESC=[155x155] Landuse (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Here is what I've tried:

  1. Convert the "Pressure" subdataset
gdal_translate -of GTiff NETCDF:"Flood.nc":Pressure test_pressure.tif
  1. Change CRS to EPSG:4326
gdal_edit.py -a_srs EPSG:4326 test_pressure.tif

I'm new to working with this kind of data so I'm wondering if there is a better/more efficient way of doing this?

If it helps here is the file as an xarray.Dataset:

<xarray.Dataset>
Dimensions:    (NX: 155, NY: 155)
Dimensions without coordinates: NX, NY
Data variables:
    X          (NX) float32 ...
    Y          (NY) float32 ...
    Pressure   (NY, NX) float32 ...
    Radar      (NY, NX) float32 ...
    Rain6h     (NY, NX) float32 ...
    RainTot    (NY, NX) float32 ...
    RH         (NY, NX) float32 ...
    U          (NY, NX) float32 ...
    V          (NY, NX) float32 ...
    Latitude   (NY, NX) float32 ...
    Longitude  (NY, NX) float32 ...
    Terrain    (NY, NX) float32 ...
    Landuse    (NY, NX) float32 ...

UPDATE:

Here is some example data. I used the xarray package to convert the Dataset into a dataframe.

         X         Y     Pressure  Radar  Rain6h  RainTot  ...         U         V   Latitude   Longitude  Terrain  Landuse
NX NY                                                      ...
0  0   0.0       0.0  1015.783508  -20.0     0.0      0.0  ...  2.894194 -0.131631  33.019382 -133.426315      0.0     17.0
   1   0.0    9000.0  1015.769226  -20.0     0.0      0.0  ...  2.387919 -0.092172  33.099922 -133.433914      0.0     17.0
   2   0.0   18000.0  1015.757568  -20.0     0.0      0.0  ...  1.781186 -0.000224  33.180466 -133.441544      0.0     17.0
   3   0.0   27000.0  1015.749756  -20.0     0.0      0.0  ...  1.024312  0.068302  33.261024 -133.449173      0.0     17.0
   4   0.0   36000.0  1015.748901  -20.0     0.0      0.0  ...  0.201696  0.096872  33.341579 -133.456818      0.0     17.0
   5   0.0   45000.0  1015.753967  -20.0     0.0      0.0  ... -0.603568  0.075843  33.422153 -133.464493      0.0     17.0

UPDATE 2: rioxarray attempt and results. The input NetCDF file can be found here.

xds is an xarray.Dataset object

I tried to use the to_raster function and got an error:

xds.rio.set_crs("epsg:4326", inplace=True)
xds["Pressure"].rio.to_raster(test.tiff)

Error: DimensionError: x dimension not found. 'set_spatial_dims()' can address this. Data variable: pressure

So I tried that next. I used multiple variables as the 'x' and 'y' parameters to set_spatial_dims() but the NY and NX combo were the only ones that didnt error:

xds.rio.set_crs("epsg:4326", inplace=True)
xds.rio.set_spatial_dims("NY", "NX", inplace=True)
xds["Pressure"].rio.to_raster(test.tiff)

This still gave me the same error as above.

I believe this might have something to do with the dimensions not having coordinates:

Dimensions:    (NX: 155, NY: 155)
Dimensions without coordinates: NX, NY
Improve this question
8
  • 2
    Is the result correct? You do not need to run gdal_edit.py command separately, you can add -a_srs EPSG:4326 into your gdal_translate command.
    – user30184
    Commented Jan 7, 2021 at 19:05
  • Thanks for the tip! I loaded the result GeoTiff file into QGIS and it looked correct, although it was black and white. Should I worry about that warning Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute, before processing? Also the GDAL docs has this note about the NetCDF driver: "This driver is intended only for importing remote sensing and geospatial datasets in form of raster images. If you want explore all data contained in NetCDF file you should use another tools.". Should I use another method to avoid any issues with using GDAL?
    – Sogun
    Commented Jan 7, 2021 at 19:15
  • Result may be OK if the georeferencing is correct and if the only band that you have has pixel values which represent the pressure correctly. What I know about NetCDF is mostly that it feels odd.
    – user30184
    Commented Jan 7, 2021 at 19:30
  • I think each column is a band. I can create a separate question to try and understand this data but it looks like there are other types of values. I believe U and V are wind arrows. I'm trying to visualize this in the most appropriate way.
    – Sogun
    Commented Jan 7, 2021 at 21:35
  • I tried opening this in QGIS and applying OpenStreetMap tiles and it looks like the conversion to EPSG:4326 didnt work. Are there any other methods I could try?
    – Sogun
    Commented Jan 8, 2021 at 0:18

2 Answers 2

Reset to default
2

This is a work in progress and is untested, but here is what you need to do:

  1. Rename the dimensions/coordinates of x/y so they are the same.
xds = xds.rename_dims({"NX": "x", "NY": "y"}).rename_vars({"X": "x", "Y": "y"}).set_coords(["x", "y"])
  1. Write the CRS to the file.
xds.rio.write_crs("<WKT or EPSG Code>", inplace=True)
  1. Write to Tiff
xds["Pressure"].rio.to_raster("test.tif")
Improve this answer
10
  • I got this error: ValueError: Cannot rename NX to longitude because longitude already exists. Try using swap_dims instead. Looks like the dimensions have to be named x and y so I tried renaming to this: .rename_dims({"NX": "x", "NY": "y"}) Then after running xds["Pressure"].rio.to_raster("test.tif") I got this error: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific. **kwargs)
    – Sogun
    Commented Jan 8, 2021 at 21:14
  • The dim and the coords for x,y need to have the same name. I made the rename dims part first to see if that would help.
    – snowman2
    Commented Jan 8, 2021 at 21:27
  • I get this now: ValuError: IndexVariable objecs must be 1-dimensional.
    – Sogun
    Commented Jan 8, 2021 at 22:00
  • Ah, that means that your data is projected and the X and Y coordinates are the ones to use. It also means you need to figure out what the projection is for your dataset.
    – snowman2
    Commented Jan 8, 2021 at 23:21
  • I'm using data that was created by someone else so I'm not sure what the projection is. If I did figure it out, what would I need to do so this conversion to GeoTIFF can work?
    – Sogun
    Commented Jan 8, 2021 at 23:43
0

I check your data and open it using Panoply, and located at west coast. If you know the bounding box <top_left_lon> <top_left_lat> <bottom_right_lon> <bottom_right_lat>, then will much easier to convert using gdal_tarnslate.

gdal_translate -of GTiff -a_ullr <top_left_lon> <top_left_lat> <bottom_right_lon> <bottom_right_lat> -a_srs EPSG:4326 input.nc output.tif

enter image description here

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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