[ncl-talk] Masking
Mary Haley
haley at ucar.edu
Tue Nov 14 08:54:46 MST 2017
I've cleaned up the script a little and added it to our shapefiles examples
page as "shapefiles_23.ncl".
See:
http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex23
--Mary
On Mon, Nov 13, 2017 at 7:30 PM, Muhammad Omer Mughal <
m.mughal1 at postgrad.curtin.edu.au> wrote:
> Hi Mary
>
>
> Many thanks for the reply. It clarifies everything.
>
>
> I needed to pass the tc2 through this array to one of my own code to
> evaluate masked UHI intensity (Temperature difference between urban and
> rural areas in a WRF output file).
>
>
> I will look into more detail in the links you sent and I appreciate your
> comments and suggestions.
>
>
>
> Muhammad Omer Mughal
> MSc BSc Mechanical Engineering
> PhD Research Scholar
> Remote Sensing and Satellite Research Group
> Department of Imaging and Applied Physics
> Curtin University
>
> Curtin University
> Tel | +61 8 9266 7962
> Fax | +61 8 9266 2377
> Mobile | 0470 237 525
>
> Email | m.mughal1 at postgrad.curtin.edu.au <m.lynch at curtin.edu.au>
> Web | http://curtin.edu.au
>
> Curtin University is a trademark of Curtin University of Technology.
> CRICOS Provider Code 00301J (WA), 02637B (NSW)
>
>
> ------------------------------
> *From:* Mary Haley <haley at ucar.edu>
> *Sent:* Tuesday, 14 November 2017 5:26:17 AM
>
> *To:* Muhammad Omer Mughal
> *Cc:* ncl-talk at ucar.edu
> *Subject:* Re: [ncl-talk] Masking
>
> Hi Muhammad,
>
> To answer your first question about what the "mask" lines are doing:
>
> The way a mask array works is that you have two 2D arrays of the exact
> same size, where one array is an array of some values, say temperature, and
> the other is an array of 0s and 1s. We call this second array a "mask
> array".
>
> When you call the "mask" function like this, using the data array and the
> mask array:
>
> tc2_mask = mask(tc2,marray,1)
>
> you are effectively telling NCL, "At the locations where marray is equal
> to 1, return the corresponding value in tc2, and otherwise return a missing
> value."
>
> What you get back, then, is an array of the same size as tc2, but with
> some of its values set to missing ( tc2 at _FillValue ).
>
> The best way to understand how mask works is to try it out on a small
> "dummy" array.
>
> Here's a script that creates a simple 4 x 3 dummy array called "t2", and a
> mask array called "marray". It masks the t2 array in two different ways,
> and prints the values in a nice format so you can see what's going on:
>
> t2 = (/(/1,2,3/),(/4,5,6/),(/7,8,9/),(/10,11,12/)/)
> t2 at _FillValue = -999
> marray = (/(/1,0,0/),(/0,1,0/),(/1,1,1/),(/0,0,1/)/)
>
> m1 = mask(t2,marray,1)
> m0 = mask(t2,marray,0)
>
> opt = True
> opt at title = "t2"
> write_matrix (t2, ("3i6"),opt)
> opt at title = "marray"
> write_matrix (marray, ("3i6"),opt)
> opt at title = "m1"
> write_matrix (m1, ("3i6"),opt)
> opt at title = "m0"
> write_matrix (m0, ("3i6"),opt)
>
> Here's the output:
>
> t2
>
> 1 2 3
> 4 5 6
> 7 8 9
> 10 11 12
>
> marray
>
> 1 0 0
> 0 1 0
> 1 1 1
> 0 0 1
>
> m1
>
> 1 -999 -999
> -999 5 -999
> 7 8 9
> -999 -999 12
>
> m0
>
> -999 2 3
> 4 -999 6
> -999 -999 -999
> 10 11 -999
>
>
>
> I'm not sure what you mean by "I need to know how to use the variables in
> other code." In the script I sent you, I showed how to mask three
> different variables on the file, using the same mask array. Is this not
> what you are talking about?
>
> I created a slightly different script showing how to create the mask,
> apply the mask to a given set of desired variables, and then how to write
> the masked variables to a new file. I then wrote a second
> "wrf_plot_orig_and_masked.ncl" script that plots the original data and
> masked data.
>
> I'm only doing this because I'm eventually going to make a new mask
> example out of all this. But again, please, I ask that you look over the
> code and try to understand it on your own. Use our functions page to look
> up any functions you don't understand:
>
> http://www.ncl.ucar.edu/Document/Functions/list_alpha.shtml
>
> Use our examples page to look at other masking and/or shapefile examples:
>
> http://www.ncl.ucar.edu/Applications/
>
> --Mary
>
>
>
> On Sat, Nov 11, 2017 at 4:59 AM, Muhammad Omer Mughal <
> m.mughal1 at postgrad.curtin.edu.au> wrote:
>
> Hi Mary
>
>
> Many thanks for the detailed email. I have worked with the code and
> produced the masked images for the variables which I need to use in my
> other code .
>
>
> However I need to know how to use the variables in other code.
>
>
> I have tried to understand the following lines but I am still struggling
> with it .Will it be possible for you to explain the following lines and how
> to output these variables to be used in another code. I would be grateful
> for your help
>
>
>
> f = addfile("wrf_mask.nc","r")
> marray = f->WRF_mask
> tc2_mask = mask(tc2,marray,1)
> lu_mask = mask(lu,marray,1)
> tc2f_mask = mask(tc2f,marray,1
>
>
>
> Muhammad Omer Mughal
> MSc BSc Mechanical Engineering
> PhD Research Scholar
> Remote Sensing and Satellite Research Group
> Department of Imaging and Applied Physics
> Curtin University
>
> Curtin University
> Tel | +61 8 9266 7962
> Fax | +61 8 9266 2377
> Mobile | 0470 237 525
>
> Email | m.mughal1 at postgrad.curtin.edu.au <m.lynch at curtin.edu.au>
> Web | http://curtin.edu.au
>
> Curtin University is a trademark of Curtin University of Technology.
> CRICOS Provider Code 00301J (WA), 02637B (NSW)
>
>
> ------------------------------
> *From:* Mary Haley <haley at ucar.edu>
> *Sent:* Friday, 10 November 2017 4:28:25 AM
>
> *To:* Muhammad Omer Mughal
> *Cc:* ncl-talk at ucar.edu
> *Subject:* Re: [ncl-talk] Masking
>
> Muhammad,
>
> A long email follows. Please read carefully, visit the links I've
> provided, and look at the two scripts I've provided.
>
> As I said in a previous email, you have curvilinear data, which means your
> data is represented by two-dimensional latitude / longitude arrays. You
> cannot use coordinate arrays with curvilinear data, as you are trying to do
> with this code:
>
> lat1d = fspan(minlat,maxlat,nlat)
> lon1d = fspan(minlon,maxlon,nlon)
> lat1d at units = "degrees_north"
> lon1d at units = "degrees_east"
> tc2a = wrf_user_getvar(a, "T2",0)
> tc2a at _FillValue = -9999
>
> tc2a!0 = "lat"
> tc2a!1 = "lon"
> tc2a&lat = lat1d
> tc2a&lon = lon1d
>
> When you use coordinate arrays via "&", this tells NCL that you have
> rectilinear data. This will NOT work for your data. To help understand the
> difference between curvilinear and rectilinear data, please see our
> "Plotting data on a map" examples page:
>
> http://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml
>
> In particular, look at these two images, where the actual lat/lon grid is
> drawn with lines:
>
> http://www.ncl.ucar.edu/Applications/Images/dataonmap_grid_1_lg.png
> http://www.ncl.ucar.edu/Applications/Images/dataonmap_grid_3_lg.png
>
> The dataonmap_grid1.ncl script shows how to plot rectilinear data, and the
> dataonma_grid_3.ncl is plotting curvilinear data.
>
> You said that the code I sent you a couple of days ago was working, so I'm
> not sure why you went back to trying to use coordinate arrays.
>
> Here's the correct way, which is just copied-and-pasted from the email I
> sent you before.
>
> nt = 0
> tc2a
>
> = wrf_user_getvar(a, "T2",
>
> nt)
> tc2a at lat2d = wrf_user_getvar(a, "lat", nt) ; This is necessary for
> shapefile masking, but do not
> tc2a at lon2d = wrf_user_getvar(a,
>
> "lon",
>
> nt) ; write these two attributes to a new NetCDF file!
> tc2a at _FillValue = -9999
>
> opt = True
> opt at debug = True
> tc2a_mask = shapefile_mask_data(tc2a,"Shapefile/malay.shp",opt)
>
>
> If you are wanting to apply the same mask to all other variables on the
> file, then what you can do is have shapefile_mask_data return a mask array,
> which is an array with 0s and 1s, where 1s represent data that falls inside
> the area, and 0s represent data that falls outside the area. You can then
> use this mask array to mask all other variables that are defined on the
> same lat/lon grid.
>
> Please see the attached wrf_gsn_5.ncl script, which I adapted from the one
> on the web. It takes a shapefile that has the Mississippi River Basin
> outline, and uses this to mask WRF data that covers the east part of the
> United States. The WRF data doesn't cover the full area of the Mississippi
> River Basin, but I wanted to use this example because I can provide you
> with the data files so you can try it yourself.
>
> The wrf_gsn_5.ncl script writes the mask array to a NetCDF file, and then
> the wrf_mask.ncl script shows how you can read this mask array back in and
> apply it to other variables on the WRF file. This is MUCH faster than
> calling shapefile_mask_data every single time. In this case,
> shapefile_mask_data takes over a minute to create the mask array. Once you
> have this array, masking data takes a fraction of a second.
>
> Try running wrf_gsn_5.ncl first, which will create several plots (to an
> X11 window). It will also write the mask array to "wrf_mask.nc". This
> script will take over a minute to run.
>
> Then, run wrf_mask.ncl, which masks the same variables, but this time it
> uses the mask array read off wrf_mask.nc to mask the data. This script
> takes a few seconds to run. It will produce a single PNG image, which I've
> attached.
>
> I've attached the four mrb.xxx files which are the files that represent
> the shapefile.
>
> To get the WRF output fie, which is larger:
>
> wget ftp://ftp.cgd.ucar.edu//incoming/wrfout_d01_2005-12-14_13:00:00
>
> or:
>
> ftp ftp.cgd.ucar.edu
> [login as anonymous]
> cd incoming
> get wrfout_d01_2005-12-14_13:00:00
> quit
>
> --Mary
>
>
>
> On Thu, Nov 9, 2017 at 6:39 AM, Muhammad Omer Mughal <
> m.mughal1 at postgrad.curtin.edu.au> wrote:
>
> Hi Marry
>
> Thank you again for the reply. I found the problem after I used CDO the
> time dimension stripped out from the data. In order to correct this problem
> I use the following
> lat=a->XLAT
> lon=a->XLONG
> minlat=min(lat)
> maxlat=max(lat)
> minlon=min(lon)
> maxlon=max(lon)
>
> nlat=129
> nlon=210
>
> lat1d = fspan(minlat,maxlat,nlat)
> lon1d = fspan(minlon,maxlon,nlon)
> lat1d at units = "degrees_north"
> lon1d at units = "degrees_east"
> tc2a = wrf_user_getvar(a, "T2",0)
> tc2a at _FillValue = -9999
>
> tc2a!0 = "lat"
> tc2a!1 = "lon"
> tc2a&lat = lat1d
> tc2a&lon = lon1d
> tc2a_mask = shapefile_mask_data(tc2a,"/data/muhdomer/NSCC/scratch/Shapef
> iles/SGP_adm0.shp",opt)
>
> Is there a way to get all variables masked using a shapefile_mask_data
> command ?
>
>
>
>
> Muhammad Omer Mughal
> MSc BSc Mechanical Engineering
> PhD Research Scholar
> Remote Sensing and Satellite Research Group
> Department of Imaging and Applied Physics
> Curtin University
>
> Curtin University
> Tel | +61 8 9266 7962
> Fax | +61 8 9266 2377
> Mobile | 0470 237 525
>
> Email | m.mughal1 at postgrad.curtin.edu.au<mailto:m.lynch at curtin.edu.au>
> Web | http://curtin.edu.au<http://curtin.edu.au/>
>
> Curtin University is a trademark of Curtin University of Technology.
> CRICOS Provider Code 00301J (WA), 02637B (NSW)
>
>
> ________________________________
> From: Mary Haley <haley at ucar.edu>
> Sent: Thursday, 9 November 2017 12:30:26 AM
> To: Muhammad Omer Mughal
> Cc: ncl-talk at ucar.edu
> Subject: Re: [ncl-talk] Masking
>
> Muhammad,
>
> Something is not right with data2.
>
> In order to correctly mask curvilinear data, your lat/lon arrays must be
> the same size as your data. If you look at the output from data2:
>
> Variable: data2
> Type: float
> Total Size: 108360 bytes
> 27090 values
> Number of Dimensions: 2
> Dimensions and sizes: [y | 129] x [x | 210]
> Coordinates:
> Number Of Attributes: 9
> lon2d : <ARRAY of 210 elements>
> lat2d : <ARRAY of 210 elements>
> stagger :
> description : TEMP at 2 M
>
> The data2 variable is 129 x 210, but the lat2d and lon2d are each 1D
> arrays of length 210. lat2d and lon2d should be 2D arrays that are the same
> size as data2.
>
> However, and this is important, attaching lat2d and lon2d attributes to
> the data variable is an NCL specific thing for plotting, regridding, and
> masking, and should NOT be used when writing new data to a file. This is
> needlessly expensive, because you likely already have the lat, lon arrays
> written as separate variables on the file.
>
> You can look at your original WRF output file to see how curvilinear data
> is generally represented on a NetCDF file. Here's an example:
>
> float T2 ( Time, south_north, west_east )
> FieldType : 104
> MemoryOrder : XY
> description : TEMP at 2 M
> units : K
> stagger :
> coordinates : XLONG XLAT XTIME
>
> Note the "coordinates" attribute attached T2. This is what tells you that
> the lat/lon arrays associated with T2 are called "XLONG" and "XLAT" on the
> file, which is what I used earlier to get the lat/lon for masking.
>
> The data2 variable longer has this coordinates attribute, so either NCO
> didn't do something correctly, or this variable was modified in some way.
>
> Please double-check that you are using NCO correctly, and don't write
> those lat2d / lon2d attributes. Also make sure that the coordinates
> attribute is included, and that you have the correct XLAT and XLONG
> variables written to the file. Once you do this, you can then mask your
> new data in the exact same way that you did for the original T2 variable.
>
> --Mary
>
>
>
>
>
> On Wed, Nov 8, 2017 at 12:35 AM, Muhammad Omer Mughal <
> m.mughal1 at postgrad.curtin.edu.au<mailto:m.mughal1 at postgrad.curtin.edu.au>>
> wrote:
> Hi Mary
>
> Thank you for the reply. The adjustment works perfectly to mask the data.
> I have another question related to this. When I use the wrfoutput file and
> use the commands lat2d I get this output
>
> Variable: data3
> Type: float
> Total Size: 108360 bytes
> 27090 values
> Number of Dimensions: 2
> Dimensions and sizes: [south_north | 129] x [west_east | 210]
> Coordinates:
> Number Of Attributes: 8
> lon2d : <ARRAY of 27090 elements>
> lat2d : <ARRAY of 27090 elements>
> coordinates : XLONG XLAT XTIME
> stagger :
> units : K
> description : TEMP at 2 M
> MemoryOrder : XY
> FieldType : 104
>
>
> When I use this to a sliced file from nco which has changed the attributes
> of the file I get the following output.
>
> Variable: data2
> Type: float
> Total Size: 108360 bytes
> 27090 values
> Number of Dimensions: 2
> Dimensions and sizes: [y | 129] x [x | 210]
> Coordinates:
> Number Of Attributes: 9
> lon2d : <ARRAY of 210 elements>
> lat2d : <ARRAY of 210 elements>
> stagger :
> description : TEMP at 2 M
> MemoryOrder : XY
> FieldType : 104
> coordinates : XLONG XLAT
> units : K
> time : 40320
>
> Also when I use the sliced file to mask the data I get the error
> (0) shapefile_mask_data:Error: not a valid rectilinear,curvilinear, or
> unstructured grid.
>
> I will appreciate your help
> Muhammad Omer Mughal
> MSc BSc Mechanical Engineering
> PhD Research Scholar
> Remote Sensing and Satellite Research Group
> Department of Imaging and Applied Physics
> Curtin University
>
> Curtin University
> Tel | +61 8 9266 7962 <+61%208%209266%207962>
> Fax | +61 8 9266 2377 <+61%208%209266%202377>
> Mobile | 0470 237 525
>
> Email | m.mughal1 at postgrad.curtin.edu.au<mailto:m.lynch at curtin.edu.au>
> Web | http://curtin.edu.au<http://curtin.edu.au/>
>
> Curtin University is a trademark of Curtin University of Technology.
> CRICOS Provider Code 00301J (WA), 02637B (NSW)
>
>
> ________________________________
> From: Mary Haley <haley at ucar.edu<mailto:haley at ucar.edu>>
> Sent: Wednesday, 8 November 2017 12:42:12 AM
> To: Muhammad Omer Mughal
> Cc: ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
> Subject: Re: [ncl-talk] Masking
>
> Muhammad,
>
> It looks like you are trying to use a rectilinear lat/lon grid for masking
> rather than the original lat/lon data associated with your "T2" variable.
>
> If your file is indeed a WRF output file, then you should use the
> XLAT/XLONG variables on the file, rather than trying to create your own
> lat/lon arrays.
>
> Your code would look something like this:
>
> nt = 0
> data
>
> = wrf_user_getvar(a, "T2",
>
> nt)
> data at lat2d = wrf_user_getvar(a, "lat", nt)
> data at lon2d = wrf_user_getvar(a,
>
> "lon",
>
> nt)
> data at _FillValue = -9999
>
> opt = True
> opt at debug = True
> data_mask = shapefile_mask_data(data,"Shapefile/malay.shp",opt)
>
>
> For an example that masks a WRF output file, see example
> shapefile_14_mask.ncl
>
> at:
>
> http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex14
>
> Note that this example is specifying "Ohio" as the specific area in the
> shapefile to mask against, but if you simply want to use the whole
> shapefile, then you do not need to set "shape_var" or "shape_names" as this
> script is doing.
>
> --Mary
>
> On Tue, Nov 7, 2017 at 6:48 AM, Muhammad Omer Mughal <
> m.mughal1 at postgrad.curtin.edu.au<mailto:m.mughal1 at postgrad.curtin.edu.au>>
> wrote:
> Hi
> I am using a WRF output file and I did the following changes "ONLY" to the
> code[ https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl] . I want
> to mask the temperature in Malaysia using the shape file malay.shp. Can
> you kindly tell me what else do I need to change as currently I am not
> able to see the data masked within Malaysia. I will be grateful for your
> help
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
>
> begin
> MASK_INSIDE = True
>
> a = addfile("./UTC_0.nc","r")
> lat=a->XLAT
> lon=a->XLONG
>
> minlat=min(lat)
> maxlat=max(lat)
> minlon=min(lon)
> maxlon=max(lon)
>
> nlat=129
> nlon=210
>
> lat1d = fspan(minlat,maxlat,nlat)
> lon1d = fspan(minlon,maxlon,nlon)
>
> data = wrf_user_getvar(a, "T2", 0)
> data at _FillValue = -9999
>
>
>
> lat1d = fspan(minlat,maxlat,nlat)
> lon1d = fspan(minlon,maxlon,nlon)
> lat1d at units = "degrees_north"
> lon1d at units = "degrees_east"
>
> ;---Attach lat/lon coordinate array information.
> data!0 = "lat"
> data!1 = "lon"
> data&lat = lat1d
> data&lon = lon1d
>
>
>
> f = addfile("Shapefile/malay.shp", "r")
> mrb_lon = f->x
> mrb_lat = f->y
> nmrb = dimsizes(mrb_lon)
>
> min_mrb_lat = min(mrb_lat)
> max_mrb_lat = max(mrb_lat)
> min_mrb_lon = min(mrb_lon)
> max_mrb_lon = max(mrb_lon)
>
>
> www.ncl.ucar.edu<https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl>
> www.ncl.ucar.edu<http://www.ncl.ucar.edu>
> ;***** ; mask_9.ncl ; ; Concepts illustrated: ; - Drawing the Mississippi
> River Basin using data from a shapefile ; - Masking a data array based on a
> geographical ...
>
>
>
> www.ncl.ucar.edu<https://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl>
> www.ncl.ucar.edu<http://www.ncl.ucar.edu>
> ;***** ; mask_9.ncl ; ; Concepts illustrated: ; - Drawing the Mississippi
> River Basin using data from a shapefile ; - Masking a data array based on a
> geographical ...
>
>
>
>
> Muhammad Omer Mughal
> MSc BSc Mechanical Engineering
> PhD Research Scholar
> Remote Sensing and Satellite Research Group
> Department of Imaging and Applied Physics
> Curtin University
>
> Curtin University
> Tel | +61 8 9266 7962
> Fax | +61 8 9266 2377
> Mobile | 0470 237 525
>
> Email | m.mughal1 at postgrad.curtin.edu.au<mailto:m.lynch at curtin.edu.au>
> Web | http://curtin.edu.au<http://curtin.edu.au/>
>
> Curtin University is a trademark of Curtin University of Technology.
> CRICOS Provider Code 00301J (WA), 02637B (NSW)
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171114/8de41f1f/attachment.html>
More information about the ncl-talk
mailing list