1

Precursor to this question: Masking NetCDF data with a shapefile in that has more than one variable Python. Note: Data files are in the previous question.

I am trying to mask my polygon but I don't know to go by it. i used the code below to plot the SIF_743 globally. But I want to now mask my area with my shapefileShapefile. How can I do that?

import xarray as xr
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
import numpy as np

# ds=xr.open_dataset("C:\proj sif\TROPOSIF_L2B_2020-01-01.nc",group= 'PRODUCT')
ds

print of the dataset output of the first code

# SIF_743 = ds['SIF_743'].values
SIF_743

#data_1d=SIF_743.ravel()
data_1d

#ds.close()
#lon = ds['longitude'].values
lon

#lat = ds['latitude'].values
lat

#dat = pd.DataFrame({'lon': lon, 'lat': lat, 'SIF_743': SIF_743})
dat

#geom = gpd.points_from_xy(dat['lon'], dat['lat'], crs=4326)
dat =gpd.GeoDataFrame(data=dat[['SIF_743']], geometry=geom)
dat

global plot of SIF_743

10
  • Please Edit the Question put your code in the body othe Question as text. Images are not legible on all devices and require anyone who wants to help to retype your code.
    – Vince
    Jan 11, 2023 at 15:44
  • @Vince please i have edited . Jan 11, 2023 at 16:07
  • It seems as if your first line of code should fail, since you did not use raw formatting or escape the backslashes. You need to put the output of the implicit print statements, so we can see the results.
    – Vince
    Jan 11, 2023 at 16:22
  • @Vince i tried to put the images of the output of the first code but it couldnt be posted. Moreover the codes that i run were all succesfull but i want to mask the a variable which is SIF_743 Jan 11, 2023 at 16:47
  • The text output could be presented as text. I repeat: Your open_dataset call should have failed.
    – Vince
    Jan 11, 2023 at 16:59

1 Answer 1

1

Here is how to extract a subset of the data in your netCDF file using geopandas:

>>> import xarray
>>> xds = xarray.open_dataset("TROPOSIF_L2B_2020-01-01.nc", group="PRODUCT")
>>> xds
<xarray.Dataset>
Dimensions:        (time: 1, n_elem: 3092691)
Coordinates:
  * time           (time) datetime64[ns] 2020-01-01
Dimensions without coordinates: n_elem
Data variables:
    delta_time     (n_elem) datetime64[ns] ...
    SIF_743        (n_elem) float32 ...
    SIF_Corr_743   (n_elem) float32 ...
    SIF_ERROR_743  (n_elem) float32 ...
    SIF_735        (n_elem) float32 ...
    SIF_Corr_735   (n_elem) float32 ...
    SIF_ERROR_735  (n_elem) float32 ...
    latitude       (n_elem) float32 ...
    longitude      (n_elem) float32 ...
>>> xds[["SIF_743", "latitude", "longitude"]]
<xarray.Dataset>
Dimensions:    (n_elem: 3092691)
Dimensions without coordinates: n_elem
Data variables:
    SIF_743    (n_elem) float32 ...
    latitude   (n_elem) float32 ...
    longitude  (n_elem) float32 ...
>>> sif_data = xds[["SIF_743", "latitude", "longitude"]]
>>> data = sif_data.to_dataframe()
>>> data
          SIF_743   latitude   longitude
n_elem                                  
0        0.762678 -85.908318   22.401075
1        0.313964 -85.953247   22.630943
2        0.170640 -86.007652   21.142496
3        0.348423 -85.998070   22.866854
4        0.201558 -86.052917   21.364904
...           ...        ...         ...
3092686  2.628601  22.160625 -159.352020
3092687  0.581590  22.146683 -159.673080
3092688  0.091324  22.153814 -159.638855
3092689  0.308466  22.160917 -159.604614
3092690  1.005781  22.167992 -159.570343

[3092691 rows x 3 columns]
>>> import geopandas
>>> geom = geopandas.points_from_xy(data.longitude, data.latitude, crs="EPSG:4326")
>>> gdf = geopandas.GeoDataFrame(data=data[['SIF_743']], geometry=geom)
>>> shapefile = geopandas.read_file("input_shapefile.shp")
>>> subset = gdf[gdf.intersects(shapefile)]
2
  • I used the subset but the data frame was empty. C:\Users\User\anaconda3\lib\site-packages\geopandas\base.py:31: UserWarning: The indices of the two GeoSeries are different.warn("The indices of the two GeoSeries are different.")C:\Users\User\AppData\Local\Temp\ipykernel_1324\2742771728.py:1: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.Use to_crs() to reproject one of the input geometries to match the CRS of the other.Left CRS: EPSG:4326.Right CRS: COMPD_CS["WGS 84+EGM96_Geoid",GEOGCS["WGS 84",DA subset = gdf[gdf.intersects(shapefile)] Jan 12, 2023 at 15:38
  • I recommend following the recommendation in the warning.
    – snowman2
    Jan 12, 2023 at 19:10

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.