Apply a shapefile to raster images in IDL

[Jenny]

Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
if you could give me any hint about the following questions:

1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?

Thank you in advance,
Jenny

[David Fanning]

Post by Jenny
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?

It's not that it is impossible to do 1 and 2 in IDL, it
certainly is possible. The problem is that it is impossible
to tell a novice IDL user how to do it. :-(

Here is a suggestion. Go to the RSI User Contribution site
and search on "shapefile". There is a shapefile viewer there
written my Jim Pendleton. I'd study that code a bit, and see if
it doesn't at least go a little ways towards answering some
of your questions.

We need a good shapefile example somewhere. Does anyone
have anything I could clean up a bit and make available
on my web page?

Cheers,

David

[Jenny]

Hi David,
Thank you very much for your reply. I've taken a quick look at Jim
Pendleton's code. you are right, as an IDL novice, I don't think I can
do it within the timeline. I guess I have to try something else, like
rasterize the shapefile in ENVI and then interpolate it to my raster
image grid...

Regards, Jenny

Post by David Fanning

Post by Jenny
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?

It's not that it is impossible to do 1 and 2 in IDL, it
certainly is possible. The problem is that it is impossible
to tell a novice IDL user how to do it. :-(
Here is a suggestion. Go to the RSI User Contribution site
and search on "shapefile". There is a shapefile viewer there
written my Jim Pendleton. I'd study that code a bit, and see if
it doesn't at least go a little ways towards answering some
of your questions.
We need a good shapefile example somewhere. Does anyone
have anything I could clean up a bit and make available
on my web page?
Cheers,
David

[Puneeth Shankar]

Hi,
while i'm able to overlay a polygon shapefile over raster(geotiff) using imap as shown below:

imap,REFORM(ortho2[0,*,*]),/order,GEOTIFF=gtiff_VI_l_shrinked,RGB_TABLE=8
imap,roi_shapefile,OVERPLOT=1

I'm unable to do the same using MAP function, which replaces the functionality of imap.My objective is to do this without using itools.Can this be done?

Regards,
Puneeth

[Kelly Dean]

It certainly can be done, but as David has pointed out, it is not
trivial. The IDL shapefile read routine requires knowledge of IDL's Objects.

I use shapefiles to add the graphics to satellite imagery using IDL's
shapefile read routines. See my image contribution to RSI's User
Contributions at www.rsinc.com. I use the county.shp to add the Colorado
county boundaries to the AVHHR imagery centered over Colorado.

To avoid the shapefile issue, I suggestion you use ArcView (or ArcGIS)
and save the shapefile as an Excel spread sheet. Then within Excel, save
the columns of latitude and longitude values as a text file and use the
IDL's OPEN/READ routines to read the text file and save the values in
separate latitude and longitude arrays. Then use the IDL's MAP_COORD
routine to convert the latitude and longitude values to window
coordinates, assuming the polar stereographic grid you desire is the
same as IDL's MAP_SET, /Stereo --- just a suggestion. However, learning
to read shapefiles directly with the IDL objects is the best and long
term why to go. The opportunities are endless.

Kelly Dean
Colorado State University
Fort Collins, Colorado

Post by Jenny
Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?
Thank you in advance,
Jenny

[Jenny]

Hi Kelly,
Thanks for your suggestion. I'd like to try it. But my polar
stereographic grid is a special one. Can I set parameters such as the
origin of the projection and the grid size (resolution) with map_set
in IDL? Sorry, I've looked the online help of map_set in IDL, but I'm
still not sure. Are P0lat and P0lon in map_set the origin of the
projection?

Many thanks, Jenny

Post by Kelly Dean
It certainly can be done, but as David has pointed out, it is not
trivial. The IDL shapefile read routine requires knowledge of IDL's Objects.
I use shapefiles to add the graphics to satellite imagery using IDL's
shapefile read routines. See my image contribution to RSI's User
Contributions at www.rsinc.com. I use the county.shp to add the Colorado
county boundaries to the AVHHR imagery centered over Colorado.
To avoid the shapefile issue, I suggestion you use ArcView (or ArcGIS)
and save the shapefile as an Excel spread sheet. Then within Excel, save
the columns of latitude and longitude values as a text file and use the
IDL's OPEN/READ routines to read the text file and save the values in
separate latitude and longitude arrays. Then use the IDL's MAP_COORD
routine to convert the latitude and longitude values to window
coordinates, assuming the polar stereographic grid you desire is the
same as IDL's MAP_SET, /Stereo --- just a suggestion. However, learning
to read shapefiles directly with the IDL objects is the best and long
term why to go. The opportunities are endless.
Kelly Dean
Colorado State University
Fort Collins, Colorado

Post by Jenny
Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?
Thank you in advance,
Jenny

[Chris Jengo]

Jenny,
This suggestion is based on three assumptions:
1. You have ENVI+IDL (sounds like it)
2. There is only one shapefile for all the files (also sounds like it)
3. Your image files are georeferenced or the lat/lon files give a
lat/lon point for every pixel in the image (known as an Internal
Geometry Model (IGM) file), or you can otherwise georeference them.

While not trivial (especially for a novice), there is a way to do this
using ENVI routines.

First, convert your shapefile to an ENVI vector file (EVF) by clicking
Vector > Open Vector File > Shapefile (I don't think there is a
documented way to do this programmatically).

Next, write an IDL program that iterates through your files and does
the following:
-Open your image and IGM files (ENVI_OPEN_FILE)
-Create a Geographic Lookup Table (GLT) file (ENVI_GLT_DOIT)
-Apply the GLT to create a georeferenced image
(ENVI_GEOREF_FROM_GLT_DOIT)
-Open the vector file (ENVI_EVF_OPEN)
-Loop through the records and extract the vector coordinates
(ENVI_EVF_READ_RECORD)
-Close the vector file (ENVI_EVF_CLOSE)
-Convert these lat/lon points to file coordinates referenced to your
georeferenced image (ENVI_CONVERT_FILE_COORDINATES)
-Create a region of interest (ENVI_CREATE_ROI)
-Populate the ROI with your points (ENVI_DEFINE_ROI)
-Get the address of every pixel within the ROI (ENVI_GET_ROI)
-Create an array with the same dimensions as your georeferenced image
(BYTARR)
-Set pixels within ROIs equal to 1 (ex: myImage[roi_addr] = 1)
-Save the mask to a temp file or memory (ENVI_WRITE_ENVI_FILE
ENVI_ENTER_DATA)
-Apply the mask to your image (ENVI_MASK_APPLY_DOIT)
-Be sure to clean up your files (ENVI_FILE_MNG)

And that's it! Nothing to it, eh? ;-) Sorry if I missed a step or
two, I was going off the top of my head. Hope this at least gets you
started...

Chris

Post by Jenny
Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?
Thank you in advance,
Jenny

[Lorenzo Busetto]

Hi Jenny,

If your images are all of the same size I think that you can solve
your problem simply by using ENVI.

1: Open one of your images and load it in a display
2: Convert your shapefile to an ENVI vector file (EVF) by clicking

Post by Chris Jengo
Vector > Open Vector File > Shapefile

3: Load your vector on the display where you have loaded your
image (In the "available vectors list" select the vector, then click
on "load" and then on "Diplay# (your_display)"
4:In the "Vector Parameters" window that appears, select "File >
Export active layer to ROI"
5: Select "basic tools > masking > build mask" and select your
display
6: Select "Options > Import ROIS", then select the ROI you created
from your vector and click on "OK". Next, choose a name for the mask
file that you want to create and click on "Apply". Now, in the
"available bands list you should see a new image, which has value 0
outside the vector and 1 inside it.
7: Now, you can use this "0-1" image as a mask for your image. Go
to Basic tools > apply mask. Select one of the images that you want to
mask, then click on "Select mask band" and select the mask image.
Click on "OK".

If you have many images and you don't want to repeat n-times the
last step you can use an IDL procedure that automatically opens the
images and apply the mask on each one.
I think that something like this should work: it's a simple
modification of the example program that you can find on the ENVI
User's guide.

pro multiple_mask

; Select input files (Select all your images)

files_list = dialog_pickfile(/READ,title ='Select input$
files',/multiple_files)

; Select the mask File

mask_file = dialog_pickfile(/read, title = 'Select Mask File')
envi_open_file, mask_file, r_fid=m_fid

; Count the number of files

num_files = N_elements (files_list)

for count = 0, num_files-1 do begin

; Selects the n- file. Output file name is
"input_file_name"+"-masked"

in_file = files_list [count]
out_file =files_list [count]+'-masked'
envi_open_file, in_file, r_fid=fid
if (fid eq -1 or m_fid eq -1) then return


; get some useful information and set the output filename.
envi_file_query, fid, ns=ns, nl=nl, nb=nb, bname=bname

; Set the keyword parameters
dims = [-1l, 0, ns-1, 0, nl-1]
pos = lindgen(nb)
m_pos = [0]

; Call the 'doit' to apply the mask

envi_mask_apply_doit, fid=fid, pos=pos, dims=dims,m_fid=m_fid,$
m_pos=m_pos, value=0, out_name=out_file, in_memory=0, $
r_fid=r_fid

endfor

end

I hadn't fully tested it, but it should work.

Hope this helps,

Lorenzo

[Jenny]

Hi Everyone,
Thank you very much for your reply. I've taken Lorenzo's suggestion,
and it works really well. This is a great group, especially to a
novice... Though just began to use IDL a couple of weeks ago, I found
it is a golden mine, worthy to dig hard :-(.

Best Regards to all!
Jenny

Post by Lorenzo Busetto
Hi Jenny,
If your images are all of the same size I think that you can solve
your problem simply by using ENVI.
1: Open one of your images and load it in a display
2: Convert your shapefile to an ENVI vector file (EVF) by clicking

Post by Chris Jengo
Vector > Open Vector File > Shapefile

3: Load your vector on the display where you have loaded your
image (In the "available vectors list" select the vector, then click
on "load" and then on "Diplay# (your_display)"
4:In the "Vector Parameters" window that appears, select "File >
Export active layer to ROI"
5: Select "basic tools > masking > build mask" and select your
display
6: Select "Options > Import ROIS", then select the ROI you created
from your vector and click on "OK". Next, choose a name for the mask
file that you want to create and click on "Apply". Now, in the
"available bands list you should see a new image, which has value 0
outside the vector and 1 inside it.
7: Now, you can use this "0-1" image as a mask for your image. Go
to Basic tools > apply mask. Select one of the images that you want to
mask, then click on "Select mask band" and select the mask image.
Click on "OK".
If you have many images and you don't want to repeat n-times the
last step you can use an IDL procedure that automatically opens the
images and apply the mask on each one.
I think that something like this should work: it's a simple
modification of the example program that you can find on the ENVI
User's guide.
pro multiple_mask
; Select input files (Select all your images)
files_list = dialog_pickfile(/READ,title ='Select input$
files',/multiple_files)
; Select the mask File
mask_file = dialog_pickfile(/read, title = 'Select Mask File')
envi_open_file, mask_file, r_fid=m_fid
; Count the number of files
num_files = N_elements (files_list)
for count = 0, num_files-1 do begin
; Selects the n- file. Output file name is
"input_file_name"+"-masked"
in_file = files_list [count]
out_file =files_list [count]+'-masked'
envi_open_file, in_file, r_fid=fid
if (fid eq -1 or m_fid eq -1) then return
; get some useful information and set the output filename.
envi_file_query, fid, ns=ns, nl=nl, nb=nb, bname=bname
; Set the keyword parameters
dims = [-1l, 0, ns-1, 0, nl-1]
pos = lindgen(nb)
m_pos = [0]
; Call the 'doit' to apply the mask
envi_mask_apply_doit, fid=fid, pos=pos, dims=dims,m_fid=m_fid,$
m_pos=m_pos, value=0, out_name=out_file, in_memory=0, $
r_fid=r_fid
endfor
end
I hadn't fully tested it, but it should work.
Hope this helps,
Lorenzo

[David Fanning]

Post by Jenny
Thank you very much for your reply. I've taken Lorenzo's suggestion,
and it works really well. This is a great group, especially to a
novice... Though just began to use IDL a couple of weeks ago, I found
it is a golden mine, worthy to dig hard :-(.

Oh, yes, lot's of great eccentrics here. You will like it. :-)

Cheers,

David

[Jenny]

Hi Everyone,
Thank you all very much. I've taken Chris+Lorenzo's suggestions, and
it works really well. Though just started to learn IDL a couple of
weeks ago, I found it is a golden mine, worthy to dig hard. :-(
With Fortran background, it seems not too hard for me to jump to IDL
direct graphics. Since the IDL object is more powerful, I think it is
time for me to learn object programming (really scared!). I'm
wondering if someone could give me some advices on how to start with
the object...

Best regards,
Jenny

Post by Lorenzo Busetto
Hi Jenny,
If your images are all of the same size I think that you can solve
your problem simply by using ENVI.
1: Open one of your images and load it in a display
2: Convert your shapefile to an ENVI vector file (EVF) by clicking

Post by Chris Jengo
Vector > Open Vector File > Shapefile

3: Load your vector on the display where you have loaded your
image (In the "available vectors list" select the vector, then click
on "load" and then on "Diplay# (your_display)"
4:In the "Vector Parameters" window that appears, select "File >
Export active layer to ROI"
5: Select "basic tools > masking > build mask" and select your
display
6: Select "Options > Import ROIS", then select the ROI you created
from your vector and click on "OK". Next, choose a name for the mask
file that you want to create and click on "Apply". Now, in the
"available bands list you should see a new image, which has value 0
outside the vector and 1 inside it.
7: Now, you can use this "0-1" image as a mask for your image. Go
to Basic tools > apply mask. Select one of the images that you want to
mask, then click on "Select mask band" and select the mask image.
Click on "OK".
If you have many images and you don't want to repeat n-times the
last step you can use an IDL procedure that automatically opens the
images and apply the mask on each one.
I think that something like this should work: it's a simple
modification of the example program that you can find on the ENVI
User's guide.
pro multiple_mask
; Select input files (Select all your images)
files_list = dialog_pickfile(/READ,title ='Select input$
files',/multiple_files)
; Select the mask File
mask_file = dialog_pickfile(/read, title = 'Select Mask File')
envi_open_file, mask_file, r_fid=m_fid
; Count the number of files
num_files = N_elements (files_list)
for count = 0, num_files-1 do begin
; Selects the n- file. Output file name is
"input_file_name"+"-masked"
in_file = files_list [count]
out_file =files_list [count]+'-masked'
envi_open_file, in_file, r_fid=fid
if (fid eq -1 or m_fid eq -1) then return
; get some useful information and set the output filename.
envi_file_query, fid, ns=ns, nl=nl, nb=nb, bname=bname
; Set the keyword parameters
dims = [-1l, 0, ns-1, 0, nl-1]
pos = lindgen(nb)
m_pos = [0]
; Call the 'doit' to apply the mask
envi_mask_apply_doit, fid=fid, pos=pos, dims=dims,m_fid=m_fid,$
m_pos=m_pos, value=0, out_name=out_file, in_memory=0, $
r_fid=r_fid
endfor
end
I hadn't fully tested it, but it should work.
Hope this helps,
Lorenzo

[Jonathan Greenberg]

Jenny (and others):

We've been developing a "bridge" between raster and vector formats
(designed, primarily, for optimized extraction of raster data from vector
coverages) -- if anyone wants to check it out, go to
http://starspan.casil.ucdavis.edu/ -- this code is open source, and we'd
love to get feedback on it. Right now its designed for unix boxes (macos x
people shouldn't have TOO much of a problem compiling it), and its not for
the unix novice since I don't think we have a makefile at this point, but if
you'd like to check it out, please do!

--j


On 10/2/04 2:29 PM, in article

Post by Jenny
Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?
Thank you in advance,
Jenny

[Adam Erickson]

Post by Jonathan Greenberg
We've been developing a "bridge" between raster and vector formats
(designed, primarily, for optimized extraction of raster data from vector
coverages) -- if anyone wants to check it out, go to
http://starspan.casil.ucdavis.edu/ -- this code is open source, and we'd
love to get feedback on it. Right now its designed for unix boxes (macos x
people shouldn't have TOO much of a problem compiling it), and its not for
the unix novice since I don't think we have a makefile at this point, but if
you'd like to check it out, please do!
--j
On 10/2/04 2:29 PM, in article

Post by Jenny
Dear Listers,
I'm a novice in IDL, not sure if my questions are totally improper. I
want to use an Arcview format shapefile as a mask, to a huge number of
raster images, which are in a polar stereographic grid format with
separate image data files and lat/lon files. I would really appreciate
1. How to read the shapefile into an IDL variable?
2. How to tailor the shapefile to my raster image grid?
3. If the above questions are not proper or impossible in IDL, would
you give me some suggestions to do the work?
Thank you in advance,
Jenny

I know this is an old message, but I suggest using QGIS for zonal statistics :) It's much more efficient. Read my reply here: https://groups.google.com/forum/#!topic/comp.lang.idl-pvwave/4E5kR9DxybQ

Cheers,

Adam