GEOREFERENCING TUTORIAL ======================= Getting Started --------------- **Note** As an aside, when you are following along in this tutorial, you should be able to use the mouse to cut any of the commands shown here, and paste them in your window. This will save time and reduce typing errors. Start up Mosaic, put it in the background hys298 >> Mosaic& Follow the links to the Georeferencing Tutorial. That should bring you to here. We will need to use a few ipw commands, so you need to have one terminal that is running on the machine vaca. hys298 >> dxterm& (In the new window on vaca ) hys298 >> rlogin vaca hys298 >> cd quinn ( _your_name_ ) Making a Grass Database ----------------------- You have three satellite images of your region that are in ipw format, red.ipw, grn.ipw and blu.ipw, one for each color band. You now want to compare this image to other real world data. This means that you need to know exactly where on earth our image is. We want to **georeference** the image data to some real world coordinate system. hys298 >> cd quinn ( _your_name_ ) hys298 >> pwd /proj/hys298/quinn hys298 >> ls blu.hst bw.pcx grn.pcx red.pcx rgb_lin.pcx blu.pcx grn.hst red.hst rgb_hist.pcx users.etx For this exercise, we're going to use the GIS grass. This GIS has a number of powerful capabilities, among them the ability to georeference images. First, however, we need to get our image data into the grass database. Let's start by creating a new directory for our database. Now, start grass with the command. hys298 >> grass.hys Fill out the form, to create a new location (image) and mapset (_your_name_) ************************************************************************* GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: imagery_______ (enter list for a list of locations) MAPSET: quinn_________ (or mapsets within a location) DATABASE: /proj/hys298/GrassDB______________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) **************************************************************************** It will ask you if you'd like to create a new mapset, choose yes. You have now told grass that you'd like to work in the hys298/GrassDB database, under the new mapset you have chosen. You are now working in the GRASS environment. The next step is to add your images into your mapset. The easiest way is to work with the raw image data. Remember, with ipw we could read and remove the headers easily. We want raw image files, so we use the commands. **NOTE: These commands only work on vaca, that's why you opened a window there** hys298 >> rmhdr red.ipw > red.raw hys298 >> rmhdr grn.ipw > grn.raw hys298 >> rmhdr blu.ipw > blu.raw Now, we can use these files in our grass environment. **NOTE: Now you're back to the GRASS environment (not on vaca)** Mapset in Location GRASS-GRID > ls *.raw blu.raw grn.raw red.raw We will use the grass command, i.tape.other, to read these files into the grass database. Mapset in Location GRASS-GRID > i.tape.other Enter tape device name: red.raw Please mount and load tape, then hit RETURN--> **************************************************************************** Please enter the following information TAPE IDENTIFICATION Davis Satellite image, red band___________________________________________ __________________________________________________________________________ IMAGE DESCRIPTION From HYS298 class_________________________________________________________ __________________________________________________________________________ __________________________________________________________________________ __________________________________________________________________________ __________________________________________________________________________ TITLE FOR THE EXTRACTED CELL FILES red_______________________________________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ***************************************************************************** GENERIC TAPE EXTRACTION tape layout 0___ number of tape files to be skipped 0___ number of records in the remaining files to be skipped 0___ number of bytes in each record to be skipped band files 1___ number of bands on the tape data format _ band sequential (BSQ) | mark one with an x x band interleaved (BIL) | 0_____ if you select BSQ format and all the bands are in a single file, enter the total number of records in the file. Otherwise enter 0 496___ length (in bytes) of the longest record on the tape 1___ blocking factor of data in the file AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) **************************************************************************** Please mark an x by the bands you want extracted x 1 AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ***************************************************************************** select a prefix/group for the band cell files to be created Enter 'list' for a list of existing imagery groups Enter 'list -f' for a verbose listing Hit RETURN to cancel request > sat ***************************************************************************** EXTRACT please select region of the image to extract start row: 1______ end row: 360____ start col: 1______ end col: 496____ ***************************************************************************** Move the created image (sat.1) to a layer named red Mapset in Location GRASS-GRID > g.rename rast=sat.1,red RENAME [sat.1] to [red] raster header category color history misc **Repeat the above steps for the two remaining bands.** Verify that you have created these files with the command: Mapset in Location GRASS-GRID > g.list rast - - - - - - - - - - - - - - - - - - - - - - - raster files available in mapset quinn: blu grn red - - - - - - - - - - - - - - - - - -- - - - - - Display these grass layers with the following commands. Mapset in Location GRASS-GRID > g.region rast=red Mapset in Location GRASS-GRID > d.mon start=x0 Graphics driver [x0] started Mapset in Location GRASS-GRID > d.rast red You will no doubt notice the image isn't nearly as nice as when you displayed the images in xv. This is because grass is pretty bad at display capabilities. In grass, combine three color images, using the command Mapset in Location GRASS-GRID > d.rgb red=red green=grn blue=blu out=rgb 100% CREATING SUPPORT FILES FOR rgb Georeferencing to exisiting Vectors ----------------------------------- Now we will create a file that lists the images that we want to rectify. These are the three single band image files, and the rgb image file we've created. We need to create a file of these filenames. We clump these files into a group that we name 'sat'. Mapset in Location GRASS-GRID > cd $LOCATION/group/sat Mapset in Location GRASS-GRID > ls REF Mapset in Location GRASS-GRID > echo "red quinn grn quinn blu quinn rgb quinn" > REF You are now ready to rectify this data to some existing location. But, you need to tell grass what that existing location is. Before doing that, you had better create that location. Since grass doesn't allow you to work on two different locations at the same time, you need to finish your current session, and then start another session. Mapset in Location GRASS-GRID > exit ************************************************************************ GRASS SESSION WRAPUP You have just finished working on mapset: The following RASTER maps belong to it: blu bw grn red rgb There are no VECTOR maps in this mapset There are no SITES maps in this mapset Shall the mapset be saved? y/n [y] y Do you wish to selectively remove data files? y/n [n]n ************************************************************************ Restart grass in the location Yolo. Remember to use _your_name_ as the mapset name. hys298 >> grass.hys ******************************************************************** GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: Yolo__________ (enter list for a list of locations) MAPSET: quinn_________ (or mapsets within a location) DATABASE: /proj/hys298/GrassDB______________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ******************************************************************** Welcome to GRASS 4.1 (Spring 1993) Update package 3 Geographic Resources Analysis Support System (GRASS) is a Trademark of U.S. Army Construction Engineering Research Laboratories (USACERL) New releases of GRASS are coordinated and produced by the Office of GRASS Integration (OGI) located at USACERL. This version running thru the Korn Shell (/usr/bin/ksh) Help is available with the command: g.help When ready to quit enter: exit Mapset in Location GRASS-GRID > exit GRASS SESSION WRAPUP You have just finished working on mapset: There are no RASTER maps in this mapset There are no VECTOR maps in this mapset There are no SITES maps in this mapset Shall the mapset be saved? y/n [y]y Do you wish to selectively remove data files? y/n [n]n You have now created a new mapset in the Yolo County location. This mapset carries the georeferencing information of that Location, in this case a UTM coordinate system representation of Yolo County. Now, return to the ungeoreferenced location, and continue working on you image data. Remember to use _your_name_ as the mapset name. hys298 >> grass.hys *************************************************************************** GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: imagery_______ (enter list for a list of locations) MAPSET: quinn_________ (or mapsets within a location) DATABASE: /proj/hys298/GrassDB______________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) *************************************************************************** Then run i.target to identify what mapset you'll use as a reference. Mapset in Location GRASS-GRID > i.target This program targets an imagery group to a GRASS database Enter group that needs a target Enter 'list' for a list of existing imagery groups Enter 'list -f' for a verbose listing Hit RETURN to cancel request > sat *********************************************************************** Please select the target LOCATION and MAPSET for group CURRENT LOCATION: imagery CURRENT MAPSET: quinn TARGET LOCATION: Yolo________________ TARGET MAPSET: quinn_______________ (enter list for a list of locations or mapsets within a location) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) *************************************************************************** This will set the target location for the georeferenced images. Now interactively create the trasnformation between the two coordinate systems, using the command, i.vpoints. Mapset in Location GRASS-GRID > d.mon select=x0 Use your mouse to make this window bigger, which you probably want to do for the next few commands that you will be using. Do not try and resize the window after you run the next few commands or you could really confuse your computer. Mapset in Location GRASS-GRID > i.vpoints *********************************************************************** Enter imagery group to be registered Enter 'list' for a list of existing imagery groups Enter 'list -f' for a verbose listing Hit RETURN to cancel request > sat *********************************************************************** This command enters you into an interactive session with grass. When you first enter this routine, a list of rasters is displayed that you may choose from. Choose the rgb image. This is the raster that you will use to locate tie points between your image, and the georefenced vectors for Yolo County. Along the bottom of your window are a number of options that you may select with your mouse. These are: QUIT | zoom | RASTER IMAGE | VECTOR MAPS | refresh | +color | ANALYZE Now use the zoom option in the bottom of the window to zoom on the map part of your image. Of the two options, the BOX zoom method is usually the best. The next step is to display some of the vectors that exist for Yolo county. Do this by using the VECTOR option in the bottom of the window. It displays a number of possible vector maps that can be displayed. Repeat this command for as many vectors as you'd like. The more vectors, the better, usually. Also, try and make the vector colors somewhat symbolic (i.e., streams color=blue). Once you've added in all the vector files, zoom in on the UC Davis campus. This takes some knowledge of where Davis is in Yolo county, and may take several tries to get it right. Now select a number of points in your image, and their corresponding location in the Yolo county vectors. You do this by clicking the left button of the mouse on a point in your raster, and then clicking the corresponding point on the vectors with the left mouse button. You need to get about ten points, and try and spread them out through the entire range of the image if possible. When enough points have been selected, use the ANALYZE option. This shows the correlation between the points and the error associated in the fit between the individual points. The highest errors are shown in red. If an error is conspicuously high, it may be that the source and destination points really don't conicide. Review the points to make sure this hasn't happened. Check the fit by choosing the OVERLAY option while in ANALYZE. This will overlay the vectors onto the raster image. Unfortunately, the vectors don't overlay onto the zoomed image, so this option isn't as useful as it could be. Once satisfied, quit this program. The points selected tie points are saved to the file $LOCATION/group/sat/POINTS. Grass will use these points for the rectification. ************************************************************************* Mapset in Location GRASS-GRID > cat $LOCATION/group/sat/POINTS # image target status # east north east north (1=ok) # 107.776619 -153.158375 607988.115361 4266215.539075 1 182.659128 -147.350746 608846.507463 4266224.059701 1 179.825640 -138.276119 608822.865672 4266312.716418 1 197.959965 -132.037313 609000.179104 4266360.000000 1 213.260802 -126.365672 609171.582090 4266401.373134 1 171.325175 -110.485075 608746.029851 4266608.238806 1 189.459500 -104.813433 608935.164179 4266649.611940 1 158.291128 -73.052239 608645.552239 4267027.880597 1 187.192710 -70.783582 608941.074627 4267021.970149 1 210.994012 -69.082090 609201.134328 4267021.970149 1 108.421733 -101.410448 608107.701493 4266714.626866 1 119.188989 -203.500000 608095.880597 4265603.462687 1 187.192710 -199.529851 608870.149254 4265662.567164 1 238.195500 -136.574627 609419.820896 4266265.432836 1 *************************************************************************** The image transformation can be performed with the i.rectify command. This command takes the image files in the selected group, along with the tie points file, and creates new images in the selected target location. Mapset in Location GRASS-GRID > i.rectify2 Enter the group containing files to be rectified Enter 'list' for a list of existing imagery groups Enter 'list -f' for a verbose listing Hit RETURN to cancel request > sat Please enter the order of the transformation you want: 1 *********************************************************************** Please select the file(s) you wish to rectify by naming an output file red@quinn . . . . . . . . . . . . . red___________ grn@quinn . . . . . . . . . . . . . grn___________ blu@quinn . . . . . . . . . . . . . blu___________ rgb@quinn . . . . . . . . . . . . . rgb___________ (enter list by any name to get a list of existing raster files) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ************************************************************************ Please select one of the following options 1. Use the current region in the target location 2. Determine the smallest region which covers the image > 2 (Then accept defaults and answer yes to all questions.) Your georeferenced images will now be created in your mapset in the Yolo location. View these new images, and overlay the exisiting vectors, by exiting your current grass session, and starting a new one in the location Yolo with mapset _your_name_ . Mapset in Location GRASS-GRID > exit ************************************************************************ GRASS SESSION WRAPUP You have just finished working on mapset: The following RASTER maps belong to it: blu bw grn red rgb There are no VECTOR maps in this mapset There are no SITES maps in this mapset Shall the mapset be saved? y/n [y] y Do you wish to selectively remove data files? y/n [n] n ************************************************************************ Now, start the new session for Location Yolo, mapset, _your_name_. hys298 >> grass.hys ************************************************************************ GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: Yolo__________ (enter list for a list of locations) MAPSET: quinn_________ (or mapsets within a location) DATABASE: /proj/hys298/GrassDB______________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ************************************************************************ Mapset in Location GRASS-GRID > d.mon select=x0; d.frame -e Mapset in Location GRASS-GRID > g.region rast=rgb Mapset in Location GRASS-GRID > d.erase Mapset in Location GRASS-GRID > d.rast rgb Mapset in Location GRASS-GRID > d.vect NghborRd color=black Try adding some more vectors in different colors. The g.list command can help in reminding you what data you have. Mapset in Location GRASS-GRID > g.list rast - - - - - - - - - - - - - - - - - - - - - - - - raster files available in mapset quinn: blu grn red rgb tblu tgrn tred trgb - - - - - - - - - - - - - - - - - - - - - - - - Mapset in Location GRASS-GRID > g.list vect - - - - - - - - - - - - - - - - - - - - - - - vector files available in mapset PERMANENT: ConctRd PrimRd boundary railroads NghborRd SecndRd ponds streams - - - - - - - - - - - - - - - - - - - - - - - Also, the d.zoom command can be used to zoom on a small section of the image. Mapset in Location GRASS-GRID > d.zoom Buttons: Left: Establish a corner Middle: Check coordinates Right: Accept region north: 4266248.13 south: 4265499.33 east: 608965.82 west: 607857.44 This region now saved as the current region. Note: run 'd.erase' for the new region to affect the graphics. Mapset in Location GRASS-GRID > d.erase Mapset in Location GRASS-GRID > d.rast rgb Mapset in Location GRASS-GRID > d.vect NghborRd color=black Note where the fit is good and where the images does not seem to fit. The errors could be the result of a number of problems, but is most likely due to errors in the TIGER data, or due to warping by the digital camera. Georeferencing to Field Data (NOT COMPLETED) -------------------------------------------- Now we will show how to do georeferencing with respect to field GPS measurements. The assumption here is that no good georeferencing information is available in any digital map form, and instead you must use GPS information that you generated in the field. We'll start with what we did in the previous section for the group 'sat' and make a new field 'field'. We could redo the steps in the previous section, but instead let's copy what we've done so far. Mapset in Location GRASS-GRID > cd $LOCATION/group Mapset in Location GRASS-GRID > ls sat Mapset in Location GRASS-GRID > cp -r sat field Mapset in Location GRASS-GRID > cd field Mapset in Location GRASS-GRID > cat POINTS # image target status # east north east north (1=ok) # 118.669094 -204.011397 608090.757180 4265594.882507 1 106.753623 -152.705490 607982.193211 4266215.248042 1 182.909022 -147.523075 608845.535248 4266220.417755 1 187.053533 -200.383706 608871.383812 4265667.258486 1 170.475487 -110.209689 608752.480418 4266597.806789 1 197.414812 -101.399584 609026.475196 4266675.352480 1 237.823799 -135.603521 609419.373368 4266261.775457 1 107.271687 -101.399584 608101.096606 4266711.540470 1 158.560017 -72.896302 608643.916449 4267016.553525 1 This file has all the points from the previous work that we've done. But Now, we want to use our GPS points instead. We'll use the above file as a template, but fill it in with our GPS points. Let's bring up our image, and try and locate where we took our field data. We will use the interactive program d.zoom to zoom into the campus, and then use d.what.rast to find the points we want. Mapset in Location GRASS-GRID > d.erase Mapset in Location GRASS-GRID > g.region rast=rgb Mapset in Location GRASS-GRID > d.rast rgb Mapset in Location GRASS-GRID > d.zoom Buttons: Left: Establish a corner Middle: Check coordinates Right: Accept region north: -57.5 south: -237.5 east: 270.5 west: 47.5 This region now saved as current region. Note: run 'd.erase' for the new region to affect the graphics. Mapset in Location GRASS-GRID > d.erase Mapset in Location GRASS-GRID > d.rast rgb Mapset in Location GRASS-GRID > d.what.rast Buttons Left: what's here Right: quit 181(E) -146(N) rgb in quinn (997) 197(E) -170(N) rgb in quinn (999) 213(E) -183(N) rgb in quinn (899) 218(E) -177(N) rgb in quinn (785) 216(E) -127(N) rgb in quinn (999) Buttons Left: what's here Right: quit 189(E) -106(N) rgb in quinn (999) Now we can use our favorite editor to put these points in the POINTS file of of field group. Now rectify the image using these field based points.