search

fmaussion / salem

Add geolocalised subsetting, masking, and plotting operations to xarray
http://salem.readthedocs.io
Other
184 stars 43 forks source link

Diagnostic variables #18

Open fmaussion opened 8 years ago

fmaussion commented 8 years ago

List of possible WRF diagnostics I implemented in IDL and that could be added to Salem:

pulsinger commented 5 years ago

Hello,

I'm using your module 'salem' for processing WRF model output. It's a helpful tool. One question though regarding diagnostic variables: if "Z" [m] is full model height, how do I get "Z" above ground i.e. about terrain ? That is "ZAG" ? As far as I've understood "Z" is calculated using the ETA/Pressure levels (i.e. bottom_top stag).

fmaussion commented 5 years ago

Thanks for your message. For getting ZAG it is enough to substract the model topography from Z. This what ZAG would do yes. I could implement it, but unfortunately we decided to deprecate the diagnostic variables altogether: https://github.com/fmaussion/salem/issues/140

I'd suggest you move to wrf-python for most of your wrf needs, maybe stay with salem for subsetting / plotting. It's a pitty because I think Salem's implementation of diagnostic variables was a great idea, but wrf-python simply has (much) more institutional and developer support!

pskunjeer commented 5 years ago

Can I export the subset data as netcdf ? This will be useful to use statistical functions in NCL

fmaussion commented 5 years ago

Can I export the subset data as netcdf ?

Yes. You can use xarray directly: http://xarray.pydata.org/en/stable/io.html#netcdf

pskunjeer commented 5 years ago

Can I export the subset data as netcdf ?

Yes. You can use xarray directly: http://xarray.pydata.org/en/stable/io.html#netcdf

Thanks for providing the link.

I have wriiten the following code to subset:

import salem

from salem.utils import get_demo_file

ds = salem.open_xr_dataset(get_demo_file('merge.nc')) ## WRF file

shdf = salem.read_shapefile(get_demo_file('AOI_PSK_WRF.shp')) ## AOI shape file

ds_roi = ds.salem.roi(shape=shdf) ## Region of Interest with all variables

ds_roi.to_netcdf('New_AOI.nc') ## Saving the netcdf file

t2 = ds_roi.T2.isel(Time=2) ## T2 at time step 2 from new AOI subset

t2.to_netcdf('New_AOI_t2.nc') ## Saving t2 AOI netcdf

t2.salem.quick_map() ## plot of t2

When i run the code there is warning and file is saved

main:11: SerializationWarning: saving variable THIS_IS_AN_IDEAL_RUN with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable ITIMESTEP with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable IVGTYP with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable ISLTYP with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable SAVE_TOPO_FROM_REAL with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable ISEEDARR_RAND_PERTURB with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable ISEEDARR_SPPT with floating point data as an integer dtype without any _FillValue to use for NaNs main:11: SerializationWarning: saving variable ISEEDARR_SKEBS with floating point data as an integer dtype without any _FillValue to use for NaNs

When i tried to open the overall netcdf file New_AOI.nc in IDV, Panoply or NCL, the file is not open.

when i tried to open the T2 netcdf file New_AOI_t2.nc it is open properly in IDV and Panoply. In NCl the values for TEMP at 2 M (K) : min=nan max=nan.

Can you please help. i think I am making mistake in saving the file. The saved file is also not in WRF format.

Regards

pskunjeer commented 5 years ago

Any updates

fmaussion commented 5 years ago

SerializationWarning: saving variable ISEEDARR_SKEBS with floating point data as an integer dtype without any _FillValue to use for NaNs

This is because when cropping, salem (or better: xarray) replaces integer values with floats. If you want to suppress these warnings you should look at the encoding of each variable: http://xarray.pydata.org/en/stable/io.html#writing-encoded-data

IDV, Panoply or NCL

this does not seem to be a problem related to salem. WRF files have a very specific format and when you crop them, other tools might have trouble to understand the projection correctly. I don't know what these tools want.

Murk89 commented 2 years ago

Hi, Is it possible to estimate the fraction of precipitation which falls as snow with salem?

fmaussion commented 2 years ago

Is it possible to estimate the fraction of precipitation which falls as snow with salem?

This depends on the microphysics scheme you are using. WRF can have at least two ways to do that, either with a SNOWFALL variable (which we would have to add to the AccumulatedVariable pattern here), or with a fraction of precifitation that falls as snow (I can't remember the name of the variable on the top of my head). What are the relevant variables in your dataset?

Murk89 commented 2 years ago

Well I have run WRF with the standard scheme, as in haven't made any changes to the microphysics. I have RAINC (mm),RAINSH(mm),RAINNC(mm) and the SNOW (snow water equivalent) vars. Thanks.

James-96 commented 2 years ago

Hi,

I'm trying to use salem.plevel to do the interpolation. A RuntimeError keeps appear, but the process isn't finished. So, this error shows again and again. Also, the same error appeared while using zlevel.

The salem version I'm using is 0.3.7. I have no idea where is wrong. Enviroment setting? Confilct with other package or some spcefic version?

The error I met is down below:

RuntimeError:
        An attempt has been made to start a new process before the
        current process has finished its bootstrapping phase.

        This probably means that you are not using fork to start your
        child processes and you have forgotten to use the proper idiom
        in the main module:

            if __name__ == '__main__':
                freeze_support()
                ...

        The "freeze_support()" line can be omitted if the program
        is not going to be frozen to produce an executable.
fmaussion commented 2 years ago

Hi @James-96, can you please open a separate issue, with a minimal working example so that we can reproduce your problem? Thanks!

James-96 commented 2 years ago

Hi, What variables are necessary for GEOPOTENTIAL and PRESSURE? I have so many wrfout files, so I'm trying to extract the necessary vars first.

fmaussion commented 2 years ago

@James-96 : https://github.com/fmaussion/salem/blob/0ae5885e4f466fc333a14544df851501f11102fd/salem/wrftools.py#L417

James-96 commented 2 years ago

Hi again, I've met a new error. 

wrf_cli = Dataset('./wind_climatology.nc')
wrf_cli = salem.open_wrf_dataset('./wind_climatology.nc')

Traceback (most recent call last):
  File "wind_plot.py", line 19, in <module>
    wrf_cli = salem.open_wrf_dataset('./wind_climatology.nc')
  File "C:\Users\73900\AppData\Roaming\Python\Python37\site-packages\salem\sio.py", line 975, in open_wrf_dataset
    nc.variables[vn] = wrftools.Unstaggerer(v)
  File "C:\Users\73900\AppData\Roaming\Python\Python37\site-packages\salem\wrftools.py", line 66, in __init__
    self.description = ncvar.description
  File "netCDF4\_netCDF4.pyx", line 4360, in netCDF4._netCDF4.Variable.__getattr__
  File "netCDF4\_netCDF4.pyx", line 4166, in netCDF4._netCDF4.Variable.getncattr
  File "netCDF4\_netCDF4.pyx", line 1407, in netCDF4._netCDF4._get_att
  File "netCDF4\_netCDF4.pyx", line 1887, in netCDF4._netCDF4._ensure_nc_success
AttributeError: NetCDF: Attribute not found

Process finished with exit code 1

netCDF4.Dataset works well on this nc file, but salem appears an AttributeError. This nc file is extracted from wrfout files by CDO, only contained U, V and other necessary vars.

fmaussion commented 2 years ago

Can you share the ncdump -h of the file?

James-96 commented 2 years ago

Can you share the ncdump -h of the file?

$ ncdump -h wind_climatology.nc netcdf wind_climatology { dimensions: XTIME = UNLIMITED ; // (12 currently) west_east_stag = 288 ; south_north = 287 ; west_east = 287 ; south_north_stag = 288 ; bottom_top = 34 ; variables: float XTIME(XTIME) ; XTIME:standard_name = "time" ; XTIME:units = "minutes since 1981-01-01 00:00:00" ; XTIME:calendar = "standard" ; XTIME:axis = "T" ; float XLONG_U(south_north, west_east_stag) ; XLONG_U:standard_name = "longitude" ; XLONG_U:long_name = "longitude" ; XLONG_U:units = "degree_east" ; XLONG_U:_CoordinateAxisType = "Lon" ; float XLAT_U(south_north, west_east_stag) ; XLAT_U:standard_name = "latitude" ; XLAT_U:long_name = "latitude" ; XLAT_U:units = "degree_north" ; XLAT_U:_CoordinateAxisType = "Lat" ; float XLONG_V(south_north_stag, west_east) ; XLONG_V:standard_name = "longitude" ; XLONG_V:long_name = "longitude" ; XLONG_V:units = "degree_east" ; XLONG_V:_CoordinateAxisType = "Lon" ; float XLAT_V(south_north_stag, west_east) ; XLAT_V:standard_name = "latitude" ; XLAT_V:long_name = "latitude" ; XLAT_V:units = "degree_north" ; XLAT_V:_CoordinateAxisType = "Lat" ; float U(XTIME, bottom_top, south_north, west_east_stag) ; U:units = "m s-1" ; U:coordinates = "XLAT_U XLONG_U" ; U:FieldType = 104 ; U:MemoryOrder = "XYZ" ; U:description = "x-wind component" ; U:stagger = "X" ; float V(XTIME, bottom_top, south_north_stag, west_east) ; V:units = "m s-1" ; V:coordinates = "XLAT_V XLONG_V" ; V:FieldType = 104 ; V:MemoryOrder = "XYZ" ; V:description = "y-wind component" ; V:stagger = "Y" ;

// global attributes: :CDI = "Climate Data Interface version 1.9.8 (https://mpimet.mpg.de/cdi)" ; :Conventions = "CF-1.6" ; :history = "Tue Mar 08 05:47:49 2022: cdo ymonmean wind.nc wind_climatology.nc\n", "Tue Mar 08 05:14:32 2022: cdo mergetime wrfout_1982_wind.nc wrfout_1983_wind.nc wrfout_1984_wind.nc wrfout_1985_wind.nc wrfout_1986_wind.nc wrfout_1987_wind.nc wrfout_1988_wind.nc wrfout_1989_wind.ncwrfout_1990_wind.nc wrfout_1991_wind.nc wrfout_1992_wind.nc wrfout_1993_wind.nc wrfout_1994_wind.nc wrfout_1995_wind.nc wrfout_1996_wind.nc wrfout_1997_wind.nc wrfout_1998_wind.nc wrfout_1999_wind.nc wrfout_2000_wind.nc Wangxj/wind.nc" ; :TITLE = " OUTPUT FROM WRF V3.9.1.1 MODEL" ; :START_DATE = "1981-01-01_00:00:00" ; :SIMULATION_START_DATE = "1981-01-01_00:00:00" ; :WEST-EAST_GRID_DIMENSION = 288 ; :SOUTH-NORTH_GRID_DIMENSION = 288 ; :BOTTOM-TOP_GRID_DIMENSION = 35 ; :DX = 24000.f ; :DY = 24000.f ; :SKEBS_ON = 0 ; :SPEC_BDY_FINAL_MU = 1 ; :USE_Q_DIABATIC = 0 ; :GRIDTYPE = "C" ; :DIFF_OPT = 1 ; :KM_OPT = 4 ; :DAMP_OPT = 3 ; :DAMPCOEF = 0.2f ; :KHDIF = 0.f ; :KVDIF = 0.f ; :MP_PHYSICS = 50 ; :RA_LW_PHYSICS = 4 ; :RA_SW_PHYSICS = 4 ; :SF_SFCLAY_PHYSICS = 2 ; :SF_SURFACE_PHYSICS = 2 ; :BL_PBL_PHYSICS = 2 ; :CU_PHYSICS = 5 ; :SF_LAKE_PHYSICS = 1 ; :SURFACE_INPUT_SOURCE = 3 ; :SST_UPDATE = 1 ; :GRID_FDDA = 0 ; :GFDDA_INTERVAL_M = 0 ; :GFDDA_END_H = 0 ; :GRID_SFDDA = 0 ; :SGFDDA_INTERVAL_M = 0 ; :SGFDDA_END_H = 0 ; :HYPSOMETRIC_OPT = 2 ; :USE_THETA_M = 0 ; :GWD_OPT = 0 ; :SF_URBAN_PHYSICS = 1 ; :SF_OCEAN_PHYSICS = 0 ; :SHCU_PHYSICS = 0 ; :MFSHCONV = 0 ; :FEEDBACK = 0 ; :SMOOTH_OPTION = 0 ; :SWRAD_SCAT = 1.f ; :W_DAMPING = 1 ; :DT = 60.f ; :RADT = 30.f ; :BLDT = 0.f ; :CUDT = 5.f ; :AER_OPT = 0 ; :SWINT_OPT = 0 ; :AER_TYPE = 1 ; :AER_AOD550_OPT = 1 ; :AER_ANGEXP_OPT = 1 ; :AER_SSA_OPT = 1 ; :AER_ASY_OPT = 1 ; :AER_AOD550_VAL = 0.12f ; :AER_ANGEXP_VAL = 1.3f ; :AER_SSA_VAL = 0.85f ; :AER_ASY_VAL = 0.9f ; :MOIST_ADV_OPT = 1 ; :SCALAR_ADV_OPT = 1 ; :TKE_ADV_OPT = 1 ; :DIFF_6TH_OPT = 0 ; :DIFF_6TH_FACTOR = 0.12f ; :OBS_NUDGE_OPT = 0 ; :BUCKET_MM = -1.f ; :BUCKET_J = -1.f ; :PREC_ACC_DT = 0.f ; :ISFTCFLX = 0 ; :ISHALLOW = 0 ; :ISFFLX = 1 ; :ICLOUD = 1 ; :ICLOUD_CU = 0 ; :TRACER_PBLMIX = 1 ; :SCALAR_PBLMIX = 0 ; :YSU_TOPDOWN_PBLMIX = 0 ; :GRAV_SETTLING = 0 ; :DFI_OPT = 0 ; :SIMULATION_INITIALIZATION_TYPE = "REAL-DATA CASE" ; :WEST-EAST_PATCH_START_UNSTAG = 1 ; :WEST-EAST_PATCH_END_UNSTAG = 287 ; :WEST-EAST_PATCH_START_STAG = 1 ; :WEST-EAST_PATCH_END_STAG = 288 ; :SOUTH-NORTH_PATCH_START_UNSTAG = 1 ; :SOUTH-NORTH_PATCH_END_UNSTAG = 287 ; :SOUTH-NORTH_PATCH_START_STAG = 1 ; :SOUTH-NORTH_PATCH_END_STAG = 288 ; :BOTTOM-TOP_PATCH_START_UNSTAG = 1 ; :BOTTOM-TOP_PATCH_END_UNSTAG = 34 ; :BOTTOM-TOP_PATCH_START_STAG = 1 ; :BOTTOM-TOP_PATCH_END_STAG = 35 ; :GRID_ID = 1 ; :PARENT_ID = 0 ; :I_PARENT_START = 1 ; :J_PARENT_START = 1 ; :PARENT_GRID_RATIO = 1 ; :CEN_LAT = 30.f ; :CEN_LON = 105.f ; :TRUELAT1 = 20.f ; :TRUELAT2 = 50.f ; :MOAD_CEN_LAT = 30.f ; :STAND_LON = 105.f ; :POLE_LAT = 90.f ; :POLE_LON = 0.f ; :GMT = 0.f ; :JULYR = 1981 ; :JULDAY = 1 ; :MAP_PROJ = 1 ; :MAP_PROJ_CHAR = "Lambert Conformal" ; :MMINLU = "MODIFIED_IGBP_MODIS_NOAH" ; :NUM_LAND_CAT = 21 ; :ISWATER = 17 ; :ISLAKE = 21 ; :ISICE = 15 ; :ISURBAN = 13 ; :ISOILWATER = 14 ; :HYBRID_OPT = -1 ; :ETAC = 0.f ; :CDO = "Climate Data Operators version 1.9.8 (https://mpimet.mpg.de/cdo)" ; }

fmaussion commented 2 years ago

Salem expects to have a description attribute. Here is a regular WRF file:

float PH(Time, bottom_top_stag, south_north, west_east) ;
        PH:least_significant_digit = 2LL ;
        PH:FieldType = 104 ;
        PH:MemoryOrder = "XYZ" ;
        PH:description = "perturbation geopotential" ;
        PH:units = "m2 s-2" ;
        PH:stagger = "Z" ;
        PH:coordinates = "XLONG XLONG_U XLAT_U XLAT_V XLONG_V XLAT" ;

I suppose one could ignore it if the attr is not there. Would you mind making a PR?

James-96 commented 2 years ago

It's interesting that this error only appears in wind files. What is PR?

fmaussion commented 2 years ago

I think it appears because CDO removes this attribute for some reason. Maybe you can prevent this behavior.

Sorry yes, a PR is https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/about-pull-requests

James-96 commented 2 years ago

I'm a beginner of coding. Maybe it's a little complex for me. I will try my best to prevent it.

fmaussion commented 2 years ago

I've asked my students if someone wants to pick this up: https://github.com/fmaussion/salem/issues/219

Let's see if someone wants to do it ;-)