[ncl-talk] how can I mask ocean data from a 4D variable using 2D oro from atmos.nc - any alternatives?

Dennis Shea shea at ucar.edu
Fri Nov 20 09:04:11 MST 2015


The 'atmos.nc' contains

Variable: oro
Number of Dimensions: 2
Dimensions and sizes:   [lat | 64] x [lon | 128]

Your data set has dimensions

        x_2 = 209 ;
        y_2 = 146 ;

Two completely different resolutions ..... The 'mask' function
documentation states that the array sizes must match
    http://www.ncl.ucar.edu/Document/Functions/Built-in/mask.shtml

Obviously, they don't match ... hence the error message.
   fatal:mask: dimension sizes  of parameter 0 and parameter 1 do not match
=====

DId you look at the 'Applications' page?  You can create a reasonably
good mask for your data sets via

http://www.ncl.ucar.edu/Applications/mask.shtml
See Example 5

Good luck

On Fri, Nov 20, 2015 at 6:59 AM, Ruksana Abedin
<ruksana.abedin at gmail.com> wrote:
> Hi,
>
> I need to mask the ocean data from my model data of temperature to see
> temperatures on land only and compare that with APHRODITES later.
>
> I tried to mask with oro from atmos file given at ncl website. Problem is
> Field16 (temperature) is 4D and oro is 2D. What can be solution for this
> problem? Can I alternatively use any other mask array for this case?
>
> Please help me to solve this problem.
>
> Current output map is attached.
>
> Giving below the file details:
>
> ncdump -h wah2_sas_1985_2007.part.seasmean.temp.nc
>
> netcdf wah2_sas_1985_2007.part.seasmean.temp {
> dimensions:
>         x = 209 ;
>         y = 145 ;
>         x_2 = 209 ;
>         y_2 = 146 ;
>         z0 = 1 ;
>         z1 = 1 ;
>         z2 = 1 ;
>         time0 = UNLIMITED ; // (1 currently)
> variables:
>         float global_longitude0(y, x) ;
>                 global_longitude0:standard_name = "longitude" ;
>                 global_longitude0:long_name = "longitude" ;
>                 global_longitude0:units = "degrees_east" ;
>                 global_longitude0:_CoordinateAxisType = "Lon" ;
>         float global_latitude0(y, x) ;
>                 global_latitude0:standard_name = "latitude" ;
>                 global_latitude0:long_name = "latitude" ;
>                 global_latitude0:units = "degrees_north" ;
>                 global_latitude0:_CoordinateAxisType = "Lat" ;
>         float global_longitude1(y_2, x_2) ;
>                 global_longitude1:standard_name = "longitude" ;
>                 global_longitude1:long_name = "longitude" ;
>                 global_longitude1:units = "degrees_east" ;
>                 global_longitude1:_CoordinateAxisType = "Lon" ;
>         float global_latitude1(y_2, x_2) ;
>                 global_latitude1:standard_name = "latitude" ;
>                 global_latitude1:long_name = "latitude" ;
>                 global_latitude1:units = "degrees_north" ;
>                 global_latitude1:_CoordinateAxisType = "Lat" ;
>         float z0(z0) ;
>                 z0:standard_name = "height" ;
>                 z0:units = "m" ;
>                 z0:positive = "up" ;
>                 z0:axis = "Z" ;
>         float z1(z1) ;
>                 z1:standard_name = "height" ;
>                 z1:units = "m" ;
>                 z1:positive = "up" ;
>                 z1:axis = "Z" ;
>         float z2(z2) ;
>                 z2:standard_name = "air_pressure" ;
>                 z2:units = "hPa" ;
>                 z2:positive = "down" ;
>                 z2:axis = "Z" ;
>         double time0(time0) ;
>                 time0:standard_name = "time" ;
>                 time0:units = "days since 1959-12-01 00:00:00" ;
>                 time0:calendar = "360_day" ;
>         float field48(time0, z0, y, x) ;
>                 field48:standard_name = "eastward_wind" ;
>                 field48:long_name = "10 METRE WIND U-COMP" ;
>                 field48:units = "m s-1" ;
>                 field48:coordinates = "global_longitude0 global_latitude0" ;
>                 field48:_FillValue = -1.073742e+09f ;
>                 field48:stash_item = "225" ;
>                 field48:stash_section = "3" ;
>                 field48:field_code = "56" ;
>                 field48:cell_method = "time: mean " ;
>         float field49(time0, z0, y, x) ;
>                 field49:standard_name = "northward_wind" ;
>                 field49:long_name = "10 METRE WIND V-COMP" ;
>                 field49:units = "m s-1" ;
>                 field49:coordinates = "global_longitude0 global_latitude0" ;
>                 field49:_FillValue = -1.073742e+09f ;
>                 field49:stash_item = "226" ;
>                 field49:stash_section = "3" ;
>                 field49:field_code = "57" ;
>                 field49:cell_method = "time: mean " ;
>         float field16(time0, z0, y_2, x_2) ;
>                 field16:standard_name = "air_temperature" ;
>                 field16:long_name = "TEMPERATURE AT 1.5M" ;
>                 field16:units = "K" ;
>                 field16:coordinates = "global_longitude1 global_latitude1" ;
>                 field16:_FillValue = -1.073742e+09f ;
>                 field16:stash_item = "236" ;
>                 field16:stash_section = "3" ;
>                 field16:field_code = "16" ;
>                 field16:cell_method = "time: maximum " ;
>         float field16_1(time0, z0, y_2, x_2) ;
>                 field16_1:standard_name = "air_temperature" ;
>                 field16_1:long_name = "TEMPERATURE AT 1.5M" ;
>                 field16_1:units = "K" ;
>                 field16_1:coordinates = "global_longitude1 global_latitude1"
> ;
>                 field16_1:_FillValue = -1.073742e+09f ;
>                 field16_1:stash_item = "236" ;
>                 field16_1:stash_section = "3" ;
>                 field16_1:field_code = "16" ;
>                 field16_1:cell_method = "time: minimum " ;
>         float field90(time0, z1, y_2, x_2) ;
>                 field90:standard_name = "precipitation_flux" ;
>                 field90:long_name = "TOTAL PRECIPITATION RATE KG/M2/S" ;
>                 field90:units = "kg m-2 s-1" ;
>                 field90:coordinates = "global_longitude1 global_latitude1" ;
>                 field90:_FillValue = -1.073742e+09f ;
>                 field90:stash_item = "216" ;
>                 field90:stash_section = "5" ;
>                 field90:field_code = "90" ;
>                 field90:cell_method = "time: mean " ;
>         float field93(time0, z1, y_2, x_2) ;
>                 field93:standard_name = "snowfall_amount" ;
>                 field93:long_name = "SNOW MASS AFTER HYDROLOGY KG/M2" ;
>                 field93:units = "kg m-2" ;
>                 field93:coordinates = "global_longitude1 global_latitude1" ;
>                 field93:_FillValue = -1.073742e+09f ;
>                 field93:stash_item = "23" ;
>                 field93:stash_section = "8" ;
>                 field93:field_code = "93" ;
>                 field93:cell_method = "time: mean " ;
>         float field1532(time0, z1, y_2, x_2) ;
>                 field1532:standard_name = "surface_runoff_flux" ;
>                 field1532:long_name = "SURFACE RUNOFF RATE KG/M2/S" ;
>                 field1532:units = "kg m-2 s-1" ;
>                 field1532:coordinates = "global_longitude1 global_latitude1"
> ;
>                 field1532:_FillValue = -1.073742e+09f ;
>                 field1532:stash_item = "234" ;
>                 field1532:stash_section = "8" ;
>                 field1532:field_code = "1532" ;
>                 field1532:cell_method = "time: mean " ;
>         float field1(time0, z2, y_2, x_2) ;
>                 field1:standard_name = "geopotential_height" ;
>                 field1:long_name = "GEOPOTENTIAL HEIGHT: PRESSURE LEVELS" ;
>                 field1:units = "m" ;
>                 field1:coordinates = "global_longitude1 global_latitude1" ;
>                 field1:_FillValue = -1.073742e+09f ;
>                 field1:stash_item = "202" ;
>                 field1:stash_section = "16" ;
>                 field1:field_code = "1" ;
>                 field1:cell_method = "time: mean " ;
>
> // global attributes:
>                 :CDI = "Climate Data Interface version 1.6.0
> (http://code.zmaw.de/projects/cdi)" ;
>                 :Conventions = "CF-1.3" ;
>                 :history = "Tue Nov 10 15:58:38 2015: cdo timmean
> -selmon,6/9 wah2_sas_1985_2007.part.nc
> wah2_sas_1985_2007.part.seasmean.temp.nc\n",
>                         "Mon Nov 09 12:02:27 2015: cdo mergetime
> wah2_sas_ga.pdi.1985.nc wah2_sas_1987_p2.nc wah2_sas_ga.pdj.1991.nc
> wah2_sas_ga.pdj.1992.nc wah2_sas_1994_p2.nc wah2_sas_ga.pdj.1995.nc
> wah2_sas_ga.pdj.1996.nc wah2_sas_ga.pdj.1997.nc wah2_sas_ga.pdj.1998.nc
> wah2_sas_ga.pdk.2000.nc wah2_sas_ga.pdk.2001.nc wah2_sas_ga.pdk.2002.nc
> wah2_sas_ga.pdk.2003.nc wah2_sas_ga.pdk.2004.nc wah2_sas_ga.pdk.2005.nc
> wah2_sas_ga.pdk.2006.nc wah2_sas_1985_2007.part.nc\n",
>                         "Fri Oct 23 12:59:32 2015: cdo mergetime
> 0049ga.pdk6dec.nc 0049ga.pdk7jan.nc 0049ga.pdk7feb.nc 0049ga.pdk7mar.nc
> 0049ga.pdk7apr.nc 0049ga.pdk7may.nc 0049ga.pdk7jun.nc 0049ga.pdk7jul.nc
> 0049ga.pdk7aug.nc 0049ga.pdk7sep.nc 0049ga.pdk7oct.nc 0049ga.pdk7nov.nc
> wah2_sas_ga.pdk.2006.nc" ;
>                 :alpham = 0.5f ;
>                 :CDO = "Climate Data Operators version 1.6.0
> (http://code.zmaw.de/projects/cdo)" ;
>
>
> while,
>
> Variable: oro
> Type: float
> Total Size: 32768 bytes
>             8192 values
> Number of Dimensions: 2
> Dimensions and sizes:   [lat | 64] x [lon | 128]
> Coordinates:
>             lat: [-87.8638..87.8638]
>             lon: [ 0..357.1875]
> Number Of Attributes: 2
>   time :         791
>   long_name :   ocean (0), land (1), sea ice (2) flag
> The script I used is:
>
> ;*********************************
> ;
> ;**********************************
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"
> ;************************************************
>
> begin
>
> ; read in data
>
> f=addfile
> ("/ouce-home/students/rege1218/WAH2_SAS_Results_1985_2011/wah2_sas_1985_2007.part.seasmean.temp.nc",
> "r");
> f1=addfile("/ouce-home/students/rege1218/WAH2_SAS_Results_1985_2011/atmos.nc","r")
>
> temp_mod1 = f->field16
> temp_mod1 = temp_mod1 - 273.15 ; Convert Kelvin -> Celsius
>
> temp_mod2 = f->field16_1
> temp_mod2 = temp_mod2 - 273.15 ; Convert Kelvin -> Celsius
>
> temp_mod = (temp_mod1+temp_mod2)/2
> time = f->time0
> lat2d = f->global_latitude1  ;0
> lon2d = f->global_longitude1 ;0
> ;lat = f->latitude1
> ;lon = f->longitude1
>
> temp_mod!0 = "time"
> temp_mod!1 = "surface"
> temp_mod!2 = "lat"
> temp_mod!3 = "lon"
> temp_mod at units = "(C)"       ; Change units to reflect.
> temp_mod at lat2d = lat2d
> temp_mod at lon2d = lon2d
>
> temp_mod = temp_mod(0,0,:,:)
> oro = f1->ORO(0,:,:)
>
> printVarSummary(temp_mod)
> printVarSummary(oro)
>
> ;***********************************************************
> ; use mask function to mask out land then ocean data
> ; ocean=0, land=1, sea_ice=2
> ;***********************************************************
> land_only = temp_mod        ;trick to keep cv's and atts
> ocean_only = temp_mod
>
>
> land_only = mask(temp_mod, oro,1)     ; Return the values of temp_mod where
> oro=1
> ocean_only = mask(temp_mod, oro,0)    ; Return the values of temp_mod where
> oro=0
>
> ;***********************************************************
> ; create plots
> ;***********************************************************
>
>   wks = gsn_open_wks("eps","temp1")
>   gsn_define_colormap(wks,"hotcolr_19lev")      ; choose colormap
>   ;gsn_reverse_colormap(wks)
>   res                          = True
>   res at cnFillOn                 = True                 ; turn on color
>   res at cnLinesOn                = False             ; turn off contour lines
>   res at cnLevelSpacingF          = 0.10              ; contour interval
>   res at cnFillDrawOrder          = "PreDraw"         ; draw contours first
>   res at lbLabelStride            = 2                 ; stride on label bar
>   res at gsnSpreadColors          = True              ; use full colormap
>   res at gsnSpreadColorEnd        = -3                ; -3 don't use land color
>   res at gsnAddCyclic             = False             ; regional data
>   res at gsnLeftString            = ""
>   res at tiMainString             = ""
>   res at gsnRightString           = ""
>   ;res at mpProjection            = "LambertConformal"
>   ;res at gsnMaskLambertConformal = True
>   res at mpMaxLatF                = 60;max(lat)  ; zoom in on region
>   res at mpMinLatF                = -30;min(lat)
>   res at mpMinLonF                = 0;min(lon)
>   res at mpMaxLonF                = 130; max(lon)
>   res at mpCenterLonF = 10                          ; def is zero
>   res at mpLandFillColor          = "Transparent"
>   res at cnLevelSelectionMode= "ManualLevels"
>   res at cnMinLevelValF      =  -4.00
>   res at cnMaxLevelValF      =  48.00 ;2
>   res at cnLevelSpacingF     = 2 ;2
>   res at lbLabelAutoStride   = True
>   temp_mod at _FillValue    = 0.0
>   temp_mod at _FillValue    = -999
>
>   plot_test = gsn_csm_contour_map(wks,land_only(0,0,:,:),res)
>
>
> end
>
> error I got:
>
> fatal:mask: dimension sizes  of parameter 0 and parameter 1 do not match
> fatal:Execute: Error occurred at or near line 54 in file
> Correctred_wah2_sas_seasonalmean_temp_1985_2007_part.ncl
>
>
> Thank you,
> Ruksana
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>


More information about the ncl-talk mailing list