Hi,
The pt_in_poly function makes a variable on the same XY grid as the
incoming variable, so it's just marking whether grid points lie
inside or outside the polygon.
Here are the commands in your script,
let xgp =
{77,78,78,81,81,82,85,85,88,88,86,86,87,87,86,86,82,82,83,83,84,84,83,83,79,79,78,77,77}
let ygp =
{31,31,30,30,29,29,28,28,27,24,24,24,23,20,20,19,19,20,20,21,21,24,24,25,25,27,27,27,27}
let gp = pt_in_poly(rf[d=1],xgp, ygp)
let cne_data_2000 = if gp eq 1 then gp*rf[d=1,l=18050:18171@ave]
First, I would suggest that you make your polygon a closed polygon
by repeating the first point at the end of the lists. This ensures
that the function will work correctly. It works by drawing lines
across the polygon and seeing how many times the line crosses the
outline of the polygon. It might work fine with the unclosed
polygon, but this will make it certain, and it means that you can
draw the polygon on a map with "plot/vs/line xgp, ygp" .
let xgp =
{77,78,78,81,81,82,85,85,88,88,86,86,87,87,86,86,82,82,83,83,84,84,83,83,79,79,78,77,77,77}
let ygp =
{31,31,30,30,29,29,28,28,27,24,24,24,23,20,20,19,19,20,20,21,21,24,24,25,25,27,27,27,27,31}
Now, I assume that what you want is only the data inside the
polygon, listed with their longitude and latitude values. You'll
need to make the result of the function into a 1-D list, define
variables with the corresponding longitude and latitude data, and
then list out just the "lon, lat, data" for the points within the
polygon.
! Define 2-D arrays, with the longitude values and the
latitude
! values at each location of the grid.
! (Try a SHADE plot of these variables to see what we are
doing here)
let xpts = x[gx=rf[d=1]] + 0*y[gy=rf[d=1]]
let ypts = 0*x[gx=rf[d=1]] + y[gy=rf[d=1]]
! 1-Dimensional lists with the locations
and all of the data in the
! variable cne_data_2000
let xlist = xsequence(xpts)
let ylist = xsequence(ypts)
let cne_data_2000_list = xsequence(cne_data_2000)
! We want just the locations where cne_data_2000
is valid.
! To pick out just some data in a variable use a mask.
let mask = if cne_data_2000_list then 1
let cne_out = mask * cne_data_2000_list
let x_out = mask * xlist
let y_out = mask * ylist
! The COMPRESS functions will move the valid values of these
! to the start of the list. Then we can list just the good data
! at the start of the 1-D lists.
! Make Ferret definitions for these variables, so the listing
will have
! a useful header. You will want to set up the right
units and a
! descriptive title for the variable.
let/title=longitude/units=degrees_east lon = compressi(x_out)
let/title=latitude/units=degrees_north lat = compressi(y_out)
let/units=" "/title=" " rf_in_region = compressi(cne_out)
let npts = mask[i=@ngd]
list/i=1:`npts` lon, lat, rf_in_region
On 1/11/2015 8:22 PM, Nitin Patil
wrote:
Dear ferret users,
In the attached script I want to extract cne_data_2000 into
ascii but what is happening is when I am doing:
list cne_data_2000 it is writing all the grid points in the
domain. I require the points/values to be listed xgp and ygp
only. Anyone knows where I am making mistake.
|