CDO

Hi Juliana,

you have to download it from https://www.ncl.ucar.edu/Applications/Scripts/shapefile_utils.ncl and e.g. save it in $HOME/NCL/shapefiles/. It is not really part of the official NCL software.

-Karin

Hi Karin,

thank you for you response. Yes, I believe I did all that but still retrieving the same issue. Please see image attached. Maybe I am doing something wrong because I've never worked with ncl before.

You can't load the file shapefile_utils.ncl in the terminal, you have to run an NCL script containing the load command with 'ncl your-script.ncl' in your terminal.

If you are not familiar with NCL read the NCL User Guide https://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide.

You can find an article about masking with shapefiles in our DKRZ vis blog https://www.dkrz.de/up/services/analysis/visualization/vis-blog/convert-an-esri-shapefile-into-a-netcdf-file.

1. ------------- REPOSTING, accidentaly deleted ------------##

Hi Karin,

thank you for you response. Yes, I believe I did all that but still retrieving the same issue. Please see image attached. Maybe I am doing something wrong because I've never worked with ncl before.

ANSWER:

You can't load the file shapefile_utils.ncl in the terminal, you have to run an NCL script containing the load command with 'ncl your-script.ncl' in your terminal.

If you are not familiar with NCL read the NCL User Guide https://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide.

You can find an article about masking with shapefiles in our DKRZ vis blog https://www.dkrz.de/up/services/analysis/visualization/vis-blog/convert-an-esri-shapefile-into-a-netcdf-file.

--Karin

New Problem!

Hi Karin,

after figuring it out how the ncl runs scripts and adapting the following script to my data, I retrieve an error (see attached picture).

How can I fix it? Thank you for your support in this.

The script looks like this:

load "/mnt/lustre02/work/mh1212/Juliana/Shapefile/shapefile_utils.ncl" ;-- load additional lib

shp_name = "./npl_admbnda_adm0_nd_20201117.shp" ;-- shapefile name
mask_name = "./npl_mask.nc" ;-- output mask file name

print_shapefile_info(shp_name)

;-- open data file to get the grid to be used
f = addfile("./tas_Amon_CESM2_historical_r1i1p1f1_gn_195001-201012_May.nc", "r")

;-- read variable
v = f->tas
v@lat = f->lat
v@lon = f->lon

;-- resources for the shapefile_mask_data function
opt = True
opt@return_mask = True ;-- this forces the return of a 0s and 1s mask array

;-- create the mask based on the given shapefile; add coordinates
mask_array = shapefile_mask_data(v, shp_name, opt)
mask_array!0 = "lat"
mask_array!1 = "lon"
mask_array@coordinates = "lat lon"
;-- mask_array&lat = f->lat
;-- mask_array&lon = f->lon
;-- mask_array@comment = "Mask shapefile: npl_admbnda_adm0_nd_20201117.shp"

;-- set value 0 to missing value
mask_array = where(mask_array .eq. 1, mask_array, mask_array@_FillValue)

;-- delete, create, and write mask array to new file
system("rm f " + mask_name)
fout = addfile(mask_name,"c")
fout>mask_array = mask_array
fout->lat = f->lat
fout->lon = f->lon

------# Thank you!

Sorry, now the code in a more fashionable way:

load "/mnt/lustre02/work/mh1212/Juliana/Shapefile/shapefile_utils.ncl"  ;-- load additional lib

shp_name  = "/mnt/lustre02/work/mh1212/Juliana/Shapefile/npl_admbnda_adm0_nd_20201117.shp" ;-- shapefile name
mask_name = "/mnt/lustre02/work/mh1212/Juliana/Shapefile/npl_mask.nc"   ;-- output mask file name

print_shapefile_info(shp_name)

;-- open data file to get the grid to be used
f = addfile("/mnt/lustre02/work/mh1212/Juliana/Shapefile/tas_Amon_CESM2_historical_r1i1p1f1_gn_195001-201012_May.nc", "r")

;-- read variable
v     = f->tas
v@lat = f->lat
v@lon = f->lon

;-- resources for the shapefile_mask_data function
opt = True
opt@return_mask = True    ;-- this forces the return of a 0s and 1s mask array

;-- create the mask based on the given shapefile; add coordinates
mask_array         =  shapefile_mask_data(v, shp_name, opt)
mask_array!0       = "lat" 
mask_array!1       = "lon" 
mask_array@coordinates = "lat lon" 
;-- mask_array&lat     =  f->lat
;-- mask_array&lon     =  f->lon
;-- mask_array@comment = "Mask shapefile: npl_admbnda_adm0_nd_20201117.shp" 

;-- set value 0 to missing value
mask_array = where(mask_array .eq. 1, mask_array, mask_array@_FillValue)

;-- delete, create, and write mask array to new file
system("rm -f " + mask_name)
fout = addfile(mask_name,"c")
fout->mask_array = mask_array
fout->lat       =  f->lat
fout->lon     =  f->lon

You are almost done. The variable tas has the dimensions tas(time, lat, lon) and shapefile_mask_data needs tas(lat,lon). You have to select one timestep first to create the mask file.

Here, we select the first timestep when reading the variable from file:

load "/mnt/lustre02/work/mh1212/Juliana/Shapefile/shapefile_utils.ncl"  ;-- load additional lib

shp_name  = "/mnt/lustre02/work/mh1212/Juliana/Shapefile/npl_admbnda_adm0_nd_20201117.shp" ;-- shapefile name
mask_name = "/mnt/lustre02/work/mh1212/Juliana/Shapefile/npl_mask.nc"   ;-- output mask file name

print_shapefile_info(shp_name)

;-- open data file to get the grid to be used
f = addfile("/mnt/lustre02/work/mh1212/Juliana/Shapefile/tas_Amon_CESM2_historical_r1i1p1f1_gn_195001-201012_May.nc", "r")

;-- read variable
v     = f->tas(0,:,:)
v@lat = f->lat
v@lon = f->lon

;-- resources for the shapefile_mask_data function
opt = True
opt@return_mask = True    ;-- this forces the return of a 0s and 1s mask array

;-- create the mask based on the given shapefile; add coordinates

mask_array         =  shapefile_mask_data(v, shp_name, opt)
mask_array!0       = "lat" 
mask_array!1       = "lon" 
mask_array@coordinates = "lat lon" 

;-- set value 0 to missing value
mask_array = where(mask_array .eq. 1, mask_array, mask_array@_FillValue)

;-- delete, create, and write mask array to new file
system("rm -f " + mask_name)
fout = addfile(mask_name,"c")
fout->mask_array = mask_array
fout->lat        =  f->lat
fout->lon        =  f->lon

The generated mask file looks like this:

netcdf npl_mask {
dimensions:
    lat = 192 ;
    lon = 288 ;
variables:
    int mask_array(lat, lon) ;
        mask_array:coordinates = "lat lon" ;
        mask_array:_FillValue = -2147483647 ;
    double lat(lat) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "latitude" ;
        lat:units = "degrees_north" ;
        lat:axis = "Y" ;
        lat:bounds = "lat_bnds" ;
    double lon(lon) ;
        lon:standard_name = "longitude" ;
        lon:long_name = "longitude" ;
        lon:units = "degrees_east" ;
        lon:axis = "X" ;
        lon:bounds = "lon_bnds" ;
}