Re: smooth boundaries of a masked area

Hi Karin,

Thank you for your response.

But this is not exactly what I meant.
Lets say I want to apply an additional wind forcing to a small region in an ocean model.
Inside some box I want the full forcing, outside of this box no forcing. At the boundaries of the box, there should be no sudden decrease but the forcing should decrease with increasing distance from the boundary.

So my first step would be to mask the data with the script you provided.
In the second step, for each gridpoint with _FillValue I would determine the distance to the next boundary of the masked box.
Then I would assign a new value to the gridpoint. For example 90% of the full forcing if the distance is 1, 80% of the full forcing if the distance is 2 and so on.

This does not seem to be very efficient and I do not know exactly how to implement this. I thought there might be some common way to do this.

Regards
Sebastian

On 19.05.2014, at 12:55, Karin Meier-Fleischer <meier-fleischer@dkrz.de> wrote:

> Hi Sebastian,
>
> to mask data along a coastline you could you shapefiles. See shapefile example mask_12.ncl:
>
> http://ncl.ucar.edu/Applications/Scripts/mask_12.ncl
> http://ncl.ucar.edu/Applications/Images/mask_12_lg.png
>
> To get the data within a choosen sub-region and masking everything else with _FillValue
> you can use the ind(..) function. The following user defined function data_mask(..)
> generates such a mask data array:
>
> ;--------------------------------------------------------------------
> undef("mask_data")
> ;--------------------------------------------------------------------
> function mask_data(data,minlat,maxlat,minlon,maxlon)
> ;--------------------------------------------------------------------
> local dims,lat,lon,nlat,nlon,lat1d,lon1d,nlatlon,ind_latlon,dims
> begin
> dims = dimsizes(data) ;-- dimension of data
> lat = data&lat ;-- retrieve latitudes
> lon = data&lon ;-- retrieve longitudes
> nlat = dims(0) ;-- dimension of lat
> nlon = dims(1) ;-- dimension of lon
>
> ;-- to use the ind function we need to convert the dimesnions lat and lon
> ;-- and the data into a 1D array
> lat1d = ndtooned(conform_dims((/nlat,nlon/),data&lat,0)) ;-- convert to 1D array
> lon1d = ndtooned(conform_dims((/nlat,nlon/),data&lon,1)) ;-- convert to 1D array
> nlatlon = dimsizes(lat1d) ;-- dimension of 1D array
> data_mask_1d = new(nlatlon,integer) ;-- assign new integer array
>
> ;-- retrieve the indices of the lat/lon's inside the sub-region
> ind_latlon = ind(lat1d.ge.minlat.and.lat1d.le.maxlat.and.lon1d.ge.minlon.and.lon1d.le.maxlon)
>
> if(.not.isatt(data,"_FillValue")) then
> data@_FillValue = default_fillvalue(typeof(data)) ;-- make sure "data" has a missing value
> end if
>
> data_1d = ndtooned(data) ;-- convert data to 1D array
> data_mask = new(dimsizes(data_1d),typeof(data),data@_FillValue) ;-- create 1D masked data array
>
> data_mask_1d = 0 ;-- initialize array to 0
>
> do n=0,dimsizes(ind_latlon)-1
> nn = ind_latlon(n) ;-- get index of point we're checking
> if(lat1d(nn).lt.minlat.or.lat1d(nn).gt.maxlat.or.lon1d(nn).lt.minlon.or.lon1d(nn).gt.maxlon)
> continue
> else
> data_mask_1d(nn) = 1 ;-- its inside the sub-region
> end if
> end do
>
> if(num(data_mask_1d.eq.1).gt.0) then
> data_mask = where(data_mask_1d.eq.1,data_1d,data_1d@_FillValue) ;-- get data or _FillValue
> else
> print("WARNING: no points inside sub-region (Lat: "+minlat+","+maxlat+" Lon: "+minlon+","+maxlon+")")
> exit
> end if
>
> ;-- for plotting purposes convert the masked data back to a 2D data array
> data_mask_2d = onedtond(data_mask,(/nlat,nlon/)) ;-- convert back to a 2d data array
> data_mask_2d!0 = "lat" ;-- set named dimension lat
> data_mask_2d!1 = "lon" ;-- set named dimension lon
> data_mask_2d&lat = lat ;-- use lat
> data_mask_2d&lon = lon ;-- use lon
> data_mask_2d&lat@units = "degrees_north" ;-- set lat units
> data_mask_2d&lon@units = "degrees_east" ;-- set lon units
>
> return(data_mask_2d) ;-- return 2D data mask array
> end
>
>
> Hope this helps,
> Karin
>
>
> Am 19.05.14 10:08, schrieb Sebastian Milinski:
>> Dear all,
>>
>> I would like to mask a region in a tripolar grid. Outside this region the values should be 0. In the transition zone (e.g. 3°) I want the values to decay (linearly or exponentially) to 0.
>> My current approach is to create a weight file with values from 0 to 1 on a regular grid. I will then interpolate it to the tripolar grid and multiply it with the data.
>>
>> Ho can I create this transition zone in NCL? For straight boundaries in N/S or E/W direction I could introduce some latitude/longitude dependence.
>> But what about more complex shapes? E.g. a masked region along a coastline?
>>
>> I thought this should be a common problem when applying a regionally confined forcing in a model. I would very much appreciate any hints.
>>
>> Regards
>> Sebastian
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
> --
> Dipl. Geophys. Karin Meier-Fleischer
> Application Support, Visualization
>
> Deutsches Klimarechenzentrum GmbH E-Mail: meier-fleischer@dkrz.de
> Bundesstrasse 45a Internet: http://www.dkrz.de/
> 20146 Hamburg Phone: +49 (0)40 460094 126
> Germany Fax: +49 (0)40 460094 270
>
> Geschäftsführer: Prof. Dr. Thomas Ludwig
> Sitz der Gesellschaft: Hamburg
> Amtsgericht Hamburg HRB 39784
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon May 19 07:02:18 2014