Re: metafile manipulation question

From: Ethan Alpert (ethan AT unknown)
Date: Thu Mar 25 1999 - 10:09:41 MST


>
> Hi,
> I need to produce contour plot of rainfall for Australian area as a
> subset from global (GCM) model and I need to have data plotted only inside
> continental outlines i.e. sea points masked out. However, GCM grid is
> rather coarse (T40 i.e. 2.8 deg) and just masking out ocean values does
> not work well. So I can produce metafile and postscript files which are
> contour of data for whole domain. The question is how can I maskout the
> colour (contour data) outside of continetal outline i.e. sea points.
> If I could produce a metafiel/PS file from say 1x1 deg where the ocean is
> white and continent is transparent then I could overlay two metafiles to
> get desired result, however I do not know how to make land data
> transparent??

        I'm not sure what version of NCAR Graphics you're using or whether
you're using the NCAR Command Language or not but here's what I think you
could do in each case. I do not believe you have to do any meta-file
manipulation.

        First, if you're using the FORTRAN low level utilities interface then
you should look into the AREAS utility. By using MAPBLA (a fortran subroutine)
and ARSCAM (FORTRAN) you can write you're own subroutine in which you
can selectively fill different map areas.

If you're doing this from scratch you'll have to look at the full documentation.

See documentation:
http://ngwww.ucar.edu/ngdoc/ng4.1/supplements/areas/
http://ngwww.ucar.edu/ngdoc/ng4.1/supplements/ezmap/

Otherwise if you're using an existing fortran program that draws filled
contours it's likely that the code is already using AREAS and EZMAPA to

You want to call mapbla to put you're map into an area map
        http://ngwww.ucar.edu/ngdoc/ng4.1/supplements/ezmap/#MAPBLA

You need to find the area and group identifiers you want to mask in the
following table:

        http://ngwww.ucar.edu/ngdoc/ng4.1/supplements/ezmap/newarealist.html

You'll then need to look for these area identifiers using ARSCAM and a
subroutine you provide. Like I said if you're already doing filled contours
only some minor modifications are needed to the routine that calls the
color fill to mask areas of the map. You just need to look for australia's
id and only call fill when the contour is in within Australia.

        http://ngwww.ucar.edu/ngdoc/ng4.1/supplements/areas/#ObtainingAreasDefined

        
        It's been a while since I used EZMAPA with AREAS and CONPACK so
I may not have given a complete picture of how to do this but in general this
is the approach.

        Now If you're using NCL (The NCAR Command Language) available in
NCAR Graphics 4.1 I've enclosed a script below that demonstrates the
masking you want to do. Before I get in to that I want to address the
1x1 degree resolution you want. NCAR Graphics 4.1 contains a SPHEREPACK 3
which can regrid using spherical harmonics. NCL has routines that will
grid from Gaussian to Gaussian, from Gaussian to Fixed Offset as well as
from Fixed to Gaussian. This method of regridding is very accurate for global
data sets.

For more information see:
http://ngwww.ucar.edu/ngdoc/ng4.1/ref/ncl/functions/g2g.html
and
http://www.scd.ucar.edu/css/software/spherepack/

Here's an NCL script that graphs some made up data but demonstrates masking.
Just save below the line and run ncl as follows:

ncl < thescript.ncl

This will create a gmeta file.

Hope this helps.

-ethan alpert

--------------------------------------------------------------------------------
begin
;
; Open an NCGM file
;
        x = create "ngcm" ncgmWorkstationClass defaultapp
                "wkColorMap" : "temp1"
        end create
;
; Make up some data so I don't have to include a data file in this message
;
        data = new((/64,128/),float)
        do i = 0,63
                data(i,:) = sin(fspan(0,2*3.141592,128))
        end do
        do j = 0,127
                data(:,j) = data(:,j) * cos(fspan(0,2*3.141592,64))
        end do
;
; Create a scalar field description needed by the contour
; to provide coordinates and data
;

        gw_lat = gaus(32)
        sf = create "sf" scalarFieldClass defaultapp
                "sfDataArray" : data
                "sfYArray": gw_lat(::-1,0)
                "sfXArray": fspan(0.0,360.0,128)
        end create

;
; Create a filled contour
;
        cn = create "cn" contourPlotClass x
                "cnScalarFieldData" : sf
                "cnFillOn" : True
                "cnFillColors" : ispan(40,63,5)
        end create

;
; Create a Lambert Confromal plot of an area encompassing Australia
;
        mp = create "mp" mapPlotClass x

                "mpProjection" : "LambertConformal"
                "mpLambertParallel1F" : -25
                "mpLambertMeridianF": 135
                "mpLambertParallel2F" : -25

                "mpLimitMode" : "CORNERS"
                "mpLeftCornerLatF" : -5
                "mpLeftCornerLonF" : 110
                "mpRightCornerLatF" : -40
                "mpRightCornerLonF" : 165
;
; Turn the fill on and set up the masking information
;

                "mpFillOn" : True
                "mpAreaMaskingOn" : True
                "mpMaskAreaSpecifiers" : "Australia"

;
; Make sure filled areas are draw after the contour plot
;

                "mpFillDrawOrder" : "postdraw"
        end create

;
; overaly the contour over the map
;
        overlay(mp,cn)
;
; draw the map and contour
;
        draw(mp)
;
; All done
;
        frame(x)
end



This archive was generated by hypermail 2b29 : Wed Jun 28 2000 - 09:40:36 MDT