1

I'm using NetCDF (TROPOSIF data) which i have been able plot and I now want to mask out my AOI with a shapefile.

I imported

import rioxarray
import geopandas as gpd

I opened the dataset using

ds = rioxarray.open_rasterio('TROPOSIF_L2B_2020-01-01.nc')
ds

Then I used ds[0] to view the xarray.Dataset

shapefile = "C:\Shapefiles\wv2modin\MODIN BY WV2.shp"
countries =gpd.read_file(shapefile)
countries

fig,ax = plt.subplots(figsize=(16,10))
countries.plot(ax=ax,column="FolderPath")

Now I can view the polygon with the above code but I don't know how to create mask because am new to Python.

link for the netcdf data shapefile

1
  • @snowman2 i used the ff code <src = rasterio.open('NETCDF:"TROPOSIF_L2B_2020-01-16.nc":/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds') lon = src.read(1)> <src = rasterio.open('NETCDF:"TROPOSIF_.-16.nc":/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds') lat = src.read(1) > <src = rasterio.open('NETCDF:"TROPOSIF_L2B_2020-01-16.nc":/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/TOA_RFL') val = src.read(1)> <dat = pd.DataFrame({'lon': lon, 'lat': lat, 'value': val}) dat> <geom = gpd.points_from_xy(dat['lon'], dat['lat'], crs=4326) dat =gpd.GeoDataFrame(data=dat[['value']], geometry=geom) dat Jan 4 at 8:20

1 Answer 1

1

https://corteva.github.io/rioxarray/stable/examples/clip_geom.html#Clip-using-a-GeoDataFrame

clipped = ds[0].rio.clip(countries.geometry.values, countries.crs, drop=False)

UPDATE:

The NetCDF file does not appear to be a raster grid:

ncdump -h TROPOSIF_L2B_2020-01-01.nc 
netcdf TROPOSIF_L2B_2020-01-01 {
dimensions:
    n_elem = 3092691 ;
    num_bd_rfl = 7 ;
    ncorner = 4 ;
    time = 1 ;

// global attributes:
        :title = "TROPOSIF_L2B" ;
        :date_created = "2021-09-14 21:26:44.262511" ;

group: METADATA {

  group: ALGORITHM_SETTINGS {

    // group attributes:
            :Polynomial\ degree\ win-743\ nm = 3LL ;
            :Number\ SVs\ win-743\ nm = 4LL ;
            :Fitting\ window\ win-743\ nm\ \(nm\) = 743., 758. ;
            :Polynomial\ degree\ win-735\ nm = 3LL ;
            :Number\ SVs\ win-735\ nm = 7LL ;
            :Fitting\ window\ win-735\ nm\ \(nm\) = 735., 758. ;
            :Cloud\ fraction\ threshold = 0.8 ;
            :SZA\ threshold = 70. ;
            :VZA\ threshold = 60. ;
            :Quality\ level\ threshold = 80LL ;
            :SIF\ reference\ wavelength\ \(nm\) = 740. ;
            :Masked-out\ spectral\ channels\ for\ SIF\ retrieval\ \(\#\) = 179LL ;
            :FWHM\ of\ macro-channels\ for\ TOA\ reflectance = 3., 3., 3. ;
    } // group ALGORITHM_SETTINGS
  } // group METADATA

group: PRODUCT {
  variables:
    float time(time) ;
        time:units = "seconds since 2010-01-01 00:00:00" ;
        time:standard_name = "time" ;
        time:comment = "Reference time of the measurements. The reference time is set to yyyy-mm-ddT00:00:00 UTC, where yyyy-mm-dd is the day on which the measurements of a particular data granule start." ;
        time:long_name = "reference start time of measurement" ;
    int delta_time(n_elem) ;
        delta_time:units = "milliseconds since 2020-01-01T00:00:00 UTC" ;
        delta_time:standard_name = "delta time" ;
        delta_time:comment = "Time difference with time for each measurement" ;
        delta_time:long_name = "offset from the reference start time of measurement" ;
    float SIF_743(n_elem) ;
        SIF_743:units = "mW/m2/sr/nm" ;
        SIF_743:standard_name = "retrieved SIF@740 743-758 nm fitting window" ;
        SIF_743:long_name = "retrieved SIF@740 (743-758nm)" ;
    float SIF_Corr_743(n_elem) ;
        SIF_Corr_743:units = "mW/m2/sr/nm" ;
        SIF_Corr_743:standard_name = "daylength-corr SIF@740 743-758 nm fitting window" ;
        SIF_Corr_743:long_name = "daylength-corr SIF@740 (743-758nm)" ;
    float SIF_ERROR_743(n_elem) ;
        SIF_ERROR_743:units = "mW/m2/sr/nm" ;
        SIF_ERROR_743:standard_name = "1-sigma error 743-758 nm fitting window" ;
        SIF_ERROR_743:long_name = "1-sigma SIF retrieval error (743-758nm)" ;
    float SIF_735(n_elem) ;
        SIF_735:units = "mW/m2/sr/nm" ;
        SIF_735:standard_name = "retrieved SIF@740 735-758 nm fitting window" ;
        SIF_735:long_name = "retrieved SIF@740 (735-758nm)" ;
    float SIF_Corr_735(n_elem) ;
        SIF_Corr_735:units = "mW/m2/sr/nm" ;
        SIF_Corr_735:standard_name = "daylength-corr SIF@740 735-758 nm fitting window" ;
        SIF_Corr_735:long_name = "daylength-corr SIF@740 (735-758nm)" ;
    float SIF_ERROR_735(n_elem) ;
        SIF_ERROR_735:units = "mW/m2/sr/nm" ;
        SIF_ERROR_735:standard_name = "1-sigma error 735-758 nm fitting window" ;
        SIF_ERROR_735:long_name = "1-sigma SIF retrieval error (735-758nm)" ;
    float latitude(n_elem) ;
        latitude:standard_name = "latitude" ;
    float longitude(n_elem) ;
        longitude:standard_name = "longitude" ;

  group: SUPPORT_DATA {

    group: DETAILED_RESULTS {
      variables:
        float TOA_RFL(n_elem, num_bd_rfl) ;
            TOA_RFL:units = "-" ;
            TOA_RFL:standard_name = "TOA Reflectance (cloud frac<0.2)" ;
            TOA_RFL:long_name = "TOA Reflectance at atmospheric windows within 665-785 nm" ;
        float WVL_RFL(num_bd_rfl) ;
            WVL_RFL:units = "nm" ;
            WVL_RFL:standard_name = "WVL_RFL" ;
            WVL_RFL:long_name = "Spectral points at which TOA_RFL is calculated" ;
        float Mean_TOA_RAD_743(n_elem) ;
            Mean_TOA_RAD_743:units = "mW/m2/sr/nm" ;
            Mean_TOA_RAD_743:standard_name = "TOA Radiance" ;
            Mean_TOA_RAD_743:long_name = "Mean TOA Radiance in 743-758 nm fitting window" ;
        float Mean_TOA_RAD_735(n_elem) ;
            Mean_TOA_RAD_735:units = "mW/m2/sr/nm" ;
            Mean_TOA_RAD_735:standard_name = "TOA Radiance" ;
            Mean_TOA_RAD_735:long_name = "Mean TOA Radiance in 735-758 nm fitting window" ;
      } // group DETAILED_RESULTS

    group: GEOLOCATIONS {
      variables:
        float viewing_zenith_angle(n_elem) ;
            viewing_zenith_angle:standard_name = "viewing zenith angle" ;
        float solar_zenith_angle(n_elem) ;
            solar_zenith_angle:standard_name = "solar zenith angle" ;
        float relative_azimuth_angle(n_elem) ;
            relative_azimuth_angle:standard_name = "relative azimuth angle" ;
        float latitude_bounds(n_elem, ncorner) ;
            latitude_bounds:standard_name = "latitude_bounds" ;
            latitude_bounds:units = "degrees_north" ;
            latitude_bounds:comment = "The four latitude boundaries of each ground pixel" ;
        float longitude_bounds(n_elem, ncorner) ;
            longitude_bounds:standard_name = "longitude_bounds" ;
            longitude_bounds:units = "degrees_east" ;
            longitude_bounds:comment = "The four longitude boundaries of each ground pixel" ;
      } // group GEOLOCATIONS

    group: INPUT_DATA {
      variables:
        float cloud_fraction_L2(n_elem) ;
            cloud_fraction_L2:standard_name = "cloud_fraction" ;
        ubyte LC_MASK(n_elem) ;
            LC_MASK:units = "([ENF=1, EBF=2, DNF=3, DBF=4, MF=5, CS=6, OS=7, WS=8, S=9, G=10, PW=11, C=12, U=13, CNV=14, SI=15, B=16])" ;
            LC_MASK:standard_name = "Land Cover Map (MODIS MCD12C1 2018)" ;
            LC_MASK:long_name = "Land Cover Map" ;
      } // group INPUT_DATA
    } // group SUPPORT_DATA
  } // group PRODUCT
}

You will need to convert this data to a grid before you can clip the data.

11
  • i tried your code but encountered an error AttributeError: 'list' object has no attribute 'rio' Jan 2 at 13:52
  • Ah, just noticed that you had multiple datasets in the response. Updated.
    – snowman2
    Jan 2 at 14:14
  • snowman2 please what could possibly be that the code now gives an error of AttributeError: 'str' object has no attribute 'rio' ? Jan 2 at 14:40
  • Can you provide the full traceback of the error?
    – snowman2
    Jan 2 at 21:41
  • snowman2 AttributeError Traceback (most recent call last) ~\AppData\Local\Temp\ipykernel_18492\3312440247.py in <module> ----> 1 clipped = ds[0].rio.clip(countries.geometry.values, countries.crs, drop=False) AttributeError: 'str' object has no attribute 'rio' Jan 3 at 7:18

Your Answer

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

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