PDA

View Full Version : PV-WAVE Mapping problem



mnatcheva
07-12-2007, 09:11 AM
Hi all,
I'm trying to display and work with maps which limits (ranges) are defined
by 4 points (form irregular polygon)? In the MAP procedure the Range
property (Range? (float) [lon1, lat1, lon2, lat2] ) is defined only by the lower
left and upper right corner of the map, assuming that the other two edges
would lay on the same longitudes/latitudes. But how should I handle the case
when the 4 end points form some totally irregular form? I found some possible solutions in IDL like below but is there any way this to be done also in the PV-WAVE:
---
stereo = MAP_PROJ_INIT('Stereographic', CENTER_LONGITUDE=-105, $
CENTER_LATITUDE=90)
; The repeat at the end will serve to close the box in oplot command
longitude = [-119.023D, -134.039, -59.959, -80.750, -119.023]
latitude = [23.117D, 53.509, 45.619, 19.8057, 23.117]
uv = MAP_PROJ_FORWARD(longitude,latitude,MAP_STRUCTURE= stereo)
---
Thanks a lot in advance!

rwagner
07-12-2007, 12:10 PM
I think the answer is probably to plot a larger square area using the map range, and then convert the lat/lon?s for the desired polygon to data coordinates, and finally to use this to mask the area desired and pull out the polygon from the image.

This is a rough solution. I'll try and work up something with some code examples in the next day or two if I can...

-Ryan

mnatcheva
07-13-2007, 06:27 AM
Hi Ryan,

Thank you very much for the fast reply! It would be very nice if you
send me also some examples - this will help me to understand you even
better. I would need to plot a map over a (polar) satellite image, for which
end points is valid the following
lower left latitude != lower right latitude
upper left latitude != upper right latitude
lower left longitude != upper left longitude
lower right longitude != upper right longitude
In this way I cannot use directly the MAP procedure with Range defined by
2 points and underlying Image. I would be very gratefull for any hints how I could overcome this!

Best regards,
M Natcheva


I think the answer is probably to plot a larger square area using the map range, and then convert the lat/lon?s for the desired polygon to data coordinates, and finally to use this to mask the area desired and pull out the polygon from the image.

This is a rough solution. I'll try and work up something with some code examples in the next day or two if I can...

-Ryan

rwagner
07-26-2007, 10:13 AM
seed = 2
xpts = (360*randomu(seed,4)-180)([0,1,2,3,0])
ypts = (90*randomu(seed,4))([0,1,2,3,0])
rnge = [ -180, min(ypts), 180, 180-min(ypts) ]
data = {,group:['cil'],area:''}
window, xsize=1000, ysize=1000
map, projection=8, range=rnge, /exact, select=data, resolution=0, /nodata
map_plots, xpts, ypts, thick=2
bmap = tvrd()
bmap = bmap eq max(bmap)
mask = blob( bmap, [0,0], [0,0] )
bdry = [ min(index_conv(bmap,where(bmap)),max=temp,dimensio n=0), temp ]
brng = bdry(1,*) - bdry(0,*)
bscl = 600.0 * brng / max(brng)
map, projection=8, range=rnge, /exact, select=data, resolution=0
plots, transpose(mask), /device, color=0, psym=3
map_plots, xpts, ypts, thick=2
imag = resamp( tvrd(bdry(0),bdry(2),brng(0),brng(1)), bscl(0), bscl(1), /i )
wait, 1
window, xsize=bscl(0), ysize=bscl(1)
tv, imag

rwagner
07-26-2007, 10:14 AM
Hi, In the previous post, I've put up an example that we've come up with. This seems to work reasonably well. Please let me know if it solves your problem.

-Ryan



P.S. It was brought to my attention that there is an unintentional space in the word dimension on the line that reads:

bdry = [ min(index_conv(bmap,where(bmap)),max=temp,dimensio n=0), temp ]

For some reason I can't get it to not show up in the post.

mnatcheva
08-14-2007, 07:18 AM
Hi Ryan,

Please excuse me for the delayed reply and thank you vey much for your
example procedure. It works really fine and fast! The only problem is that
we attempted to plot the selection of coatlines over an already projected
image. Do you know if this is also possible?

Best regards and thanks in advance,
MNatcheva

rwagner
08-14-2007, 02:24 PM
Try taking a look at the image keyword in the map procedure:

Image?Specifies an image (2D array) to be warped around the map projection and displayed with the map lines. If the image array is not of type byte, it will be automatically passed through the BYTSCL procedure before being plotted. The image will be warped to exactly cover the area specified using the Range keyword.

Perhaps you can do a TVRD on the projection, then use the result as input for the image keyword in the call to map?

-Ryan