Mask NETCDF file by shapefile?

hello everyone, could someone help me?
I'm trying to make a mask for a netcdf file just from my study area, which is the shapefile area that I have.
I used an NCL script to transform my shapefile file into netcdf. So you can use it in the CDO to make the mask. However, the CDO cannot read the file I created with the NCL
the script I used was this one, proposed by a user of this forum

load "/usr/local/lib/ncarg/nclscripts/shapefile_utils.ncl" ; 
    begin
        shpname  = "/home/Luiz/bacia_ama_bolivia.shp" 
        maskname = "/home/Luiz/mask_bolivia.nc" 
        print_shapefile_info(shpname)
        ;-- open data file to get the grid to be used
        f = addfile("GPM_DJF.nc","r")
        ;print(f)
        vNames = getfilevarnames (f) ; Pegar o nome das variaveis no arquivo
        ;print(vNames) ; mostrar o nome das variaveis
        ; Devem estar os nomes de cada variavel do arquivo selecionado
        ;-- read variable
        var                    =  f->GPM_3IMERGM_06_precipitation(:,:)
        var@lat2d              =  f->lat_bnds
        var@lon2d              =  f->lon_bnds
        ;-- shapefile mask resources
        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
        mask_array             =  shapefile_mask_data(var, shpname, opt)
        mask_array!0           = "y" 
        mask_array!1           = "x" 
        mask_array@coordinates = "lat_bnds lon_bnds" 
        ;-- create new netCDF file and write mask array
        system("rm f " + maskname)
        fout = addfile(maskname,"c")
        fout>mask_array =  mask_array
        fout->lat_bnds      =  f->lat_bnds
        fout->lon_bnds     =  f->lon_bnds
        df = addfile("mask_bolivia.nc","r")
        print(df)
    end

what am i doing wrong? can anybody help me

Download all files

bacia_ama_bolivia.dbf (314 Bytes) bacia_ama_bolivia.dbf
bacia_ama_bolivia.cpg (5 Bytes) bacia_ama_bolivia.cpg
bacia_ama_bolivia.prj (151 Bytes) bacia_ama_bolivia.prj
bacia_ama_bolivia.sbn (132 Bytes) bacia_ama_bolivia.sbn
bacia_ama_bolivia.sbx (116 Bytes) bacia_ama_bolivia.sbx
bacia_ama_bolivia.shp.xml (8 KB) bacia_ama_bolivia.shp.xml
bacia_ama_bolivia.shx (108 Bytes) bacia_ama_bolivia.shx
bacia_ama_bolivia.shp (362 KB) bacia_ama_bolivia.shp	shapefile (my study area)
GPM_DJF.nc (198 KB) GPM_DJF.nc	netcdf that needs to be masked by the shapefile

Replies (3)

RE: Mask NETCDF file by shapefile? - Added by Karin Meier-Fleischer about 4 years ago

Hi Luiz,

why are you using the lat_bnds/lon_bnds? You ahve to use lat and lon as usual.

load "/usr/local/lib/ncarg/nclscripts/shapefile_utils.ncl" ; 

begin
    shpname  = "/home/Luiz/bacia_ama_bolivia.shp" 
    maskname = "/home/Luiz/mask_bolivia.nc" 

    print_shapefile_info(shpname)

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

    vNames = getfilevarnames (f) ; Pegar o nome das variaveis no arquivo

    ;-- read variable
    var                    =  f->GPM_3IMERGM_06_precipitation(:,:)
    var@lat2d              =  f->lat
    var@lon2d              =  f->lon

    ;-- shapefile mask resources
    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
    mask_array             =  shapefile_mask_data(var, shpname, opt)
    mask_array!0           = "y" 
    mask_array!1           = "x" 
    mask_array@coordinates = "lat lon" 

    ;-- create new netCDF file and write mask array
    system("rm -f " + maskname)
    fout = addfile(maskname,"c")

    fout->mask_array =  mask_array
    fout->lat_bnds   =  f->lat_bnds
    fout->lon_bnds   =  f->lon_bnds
end

-Karin

RE: Mask NETCDF file by shapefile? - Added by L. Santos about 4 years ago

Hi Karin..It worked..was using the wrong variable. thank you very much ... and thank you for your script in ncl. helped me a lot.

RE: Mask NETCDF file by shapefile? - Added by Alessandro Ceppi about 2 years ago

Hi Karin,

I have the same problem, but with netcdf from here: https://cds.climate.copernicus.eu/cdsapp#!/dataset/sis-urban-climate-cities?tab=form
I cannot sort something out.
Can you help me, please?
here it is my script:

load "$NCARG_ROOT/lib/ncarg/nclscripts/shapefile_utils.ncl" ;

begin
shpname = "/home/lucacerri/sol_img/mappe_meteo/sol/corsico/corsico.shp"
maskname = "/home/lucacerri/sol_img/mappe_meteo/sol/corsico/corsico_mask.nc"

print_shapefile_info(shpname)

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

vNames = getfilevarnames (f) ;

;-- read variable
    var                    =  f->tas(0,:,:)
    var@lat2d              =  f->latitude
    var@lon2d              =  f->longitude

;-- shapefile mask resources
    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
    mask_array             =  shapefile_mask_data(var, shpname, opt)
    mask_array!0           = "y" 
    mask_array!1           = "x" 
    mask_array@coordinates = "latitude longitude"

;-- create new netCDF file and write mask array
    system("rm -f " + maskname)
    fout = addfile(maskname,"c")

fout->mask_array =  mask_array
    fout->latitude    =  f->latitude
    fout->longitude   =  f->longitude
end

and here are the files nectdf and shapefiles.
Many thanks in advance.
Alessandro

Download all files

corsico.rar (4.95 KB) corsico.rar	shapefile
uhi.nc (500 KB) uhi.nc	netcdf

(1-3/3)

Project

General

Profile

CDO