Re: No-frills contour plot

Hi Oli,

I've done something similar to what you're after, although I use
additional command-line tools to extract just the "map" portion of a
plot and inject georeferencing information. I first find where the
"map" lies on the postscript page, and then use Imagemagick's
"convert" tool to crop the page at those postscript coordinates. I
then use GDAL's "gdal_translate" utility to inject georeferencing
info, creating a geotiff. The whole thing can be scripted from within
an NCL script by making calls out to the "system()" function to invoke
the command-line utilities with computed values.

The location of the map on a postscript page can be determined using
the following solution provided by Dave Brown. We first find the
location on the page of the normalized device coordinate space (NDC),
using the workstation's resources:

        dc = new(4,integer)
        getvalues wks
           "wkDeviceLowerX" : dc(0)
           "wkDeviceLowerY" : dc(1)
           "wkDeviceUpperX" : dc(2)
           "wkDeviceUpperY" : dc(3)
        end getvalues

The location of the "map" within NDC space is obtained from the plot's
resources:

       vp = new(4,float)
       getvalues plot
             "vpXF" : vp(0)
             "vpYF" : vp(1)
             "vpWidthF" : vp(2)
             "vpHeightF" : vp(3)
       end getvalues

 From these, its straightforward to calculate the *location on the
page* of the map portion of the plot:

           dvp = new(4,float)
               dvp(0) = dc(0) + vp(0) * (dc(2) - dc(0))
         dvp(1) = dc(1) + (vp(1) - vp(3)) * (dc(3) - dc(1))
         dvp(2) = dc(0) + (vp(0) + vp(2)) * (dc(2) - dc(0))
         dvp(3) = dc(1) + vp(1) * (dc(3) - dc(1))

These are postscript page-coordinates that can be used as parameters
to the "-crop" switch of the convert command.

I find the geographical coordinates of the map portion of a plot by:

    getvalues plot
        "vpXF": xf;
        "vpYF": yf;
        "vpWidthF": width;
        "vpHeightF": height;
    end getvalues

    slop = .0000001
    xNDC = (/ xf+slop, xf+width-slop, xf+width-slop, xf+slop /)
    yNDC = (/ yf-slop, yf-slop, yf-height+slop, yf-height+slop /)
    xWorld = new(dimsizes(xNDC), float)
    yWorld = new(dimsizes(yNDC), float)
    ndctodata(plot, xNDC, yNDC, xWorld, yWorld)

The "slop" term in the expressions above is a work-around to the fact
that ndctodata() fails if the NDC coordinates are *exactly* on the
boundary of "plot" -- this code forces the coordinates to be slightly
inside the bounds of "plot". This obviously introduces error, but
generally of negligible magnitude. xWorld/yWorld now have geographical
coordinates.

Hope this helps...

Rick Brownrigg

From: <Oliver.Fuhrer_at_meteoswiss.ch>
> Date: March 25, 2009 11:49:33 AM MDT
> To: <ncl-talk_at_ucar.edu>
> Subject: No-frills contour plot
>
> Hi all,
>
> I would like to produce an EPS file of defined size of a contour plot
> with absolutely no information (!) except the filled contours. I've
> been
> trying to play around with options (tiMainOn=False, tmYROn=False,
> tmBorderThicknessF=0.0, cnGridBoundPerimOn=False, etc...) but did not
> really manage. First of all the EPS always contained a white border
> (see
> attached image), second I don't know how to exactly control the size
> (for example in inch) of the resulting plot. The goal would be to
> create
> a simple raster image by conversion that could easily be geo-
> referenced
> and used in other software.
>
> Has anyone already done something like this?
>
> Thanks,
> Oli
>
>
> ________________________________________
>
> Oliver Fuhrer
> Numerical Models
>
> Federal Departement of Home Affairs FDHA
> Federal Office of Meteorology and Climatology MeteoSwiss
>
> Kraehbuehlstrasse 58, P.O. Box 514, CH-8044 Zurich, Switzerland
>
> Tel. +41 44 256 93 59
> Fax +41 44 256 92 78
> oliver.fuhrer_at_meteoswiss.ch
> www.meteoswiss.ch - First-hand information
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 26 2009 - 14:30:26 MDT