Converting Altas LIDAR DEMS into Friendlier Format
Quick instructions on how to convert (hopefully) the ATLAS LIDAR DEMs, encoded in a USGS format, into a friendlier (but maybe not as rich) ASCII format.
Based on usage of GRASS GIS 6.3.
Okay. I'm making some assumptions in writing this guide.
- You know how to invoke Grass GIS
- You know how to create a new Project location and MAPSET (see this file).
Also, currently, I'll assume that you are interested in just converting a single 1/4 quad. I'll discuss how to import/convert multiple 1/4 quads later in this page.
Scenario 1. Importing/Converting Single 1/4 Quad.
1) Getting the datafile.Download the 1/4 quad of interest. To obtain this, go to the http://atlas.lsu.edu web-site. Click on download Data and scroll to the bottom of the page. Select the link that says "LIDAR" data.
Okay. Now, you can select a Parish to zoom into. Do so.
Then, select the quad of interest.
A side-frame should be on the left. This will allow you to download a 1/4 quad region. You can download DEM, contour, and so forth. Download the DEM.
Now, unzip the file. The contents will appear in the directory the zip file is located. No new directories are created. You should have a dem file, a readme, and a bunch of files (all containing about the same info) describing, to an extent, the dem file.
You only need the dem file.
2) Gathering some information and importing into GRASS GISHere, we need a program call 'gdalinfo'. Realistically, you should already have this program, as the GDAL library is needed in order for GRASS GIS to work. NOTE : commands will be italic and highlighted. NOTE: This is a command-line program.
Type: gdalinfo name_of_file.dem
You should have a lot of information casacading down the screen.
The items of importance are PIXEL_SIZE and Corner Coordinates. You will need those to create the GRASS location and MAPSET. If you don't know, please go HERE to find out how. IMPORTANT: Make sure the mapset is 'large' enough spatially to include the 1/4 quad (in other words, the corner coordinates of the dem should be contained with the corner coordinates of the MAPSET).
Assuming you have a valid Mapset that meets the requirements, the next step is to import the quad. So, start Grass
grass63
A gui should open; select the location and MAPSET that you created (or are going to use). Once that is done, you will need to run this command:
r.in.gdal input=DEM_FILE output=GRASS_NAME
where
- DEM_FILE is the file to upload, including extension
- GRASS_NAME is how the file should be known within grass. I'd suggest having a meaningful extension.
Now, some notes. You should be aware of these.
- If DEM_FILE is not in the same directory as where you executed the grass63 command, you need to also include as part of DEM_FILE the relative path.
- GRASS_NAME. Grass GIS makes an actual copy of the file (and generates some internal support info). So, make sure you know what the name means. You can make it the same as the DEM_FILE (minus any path information).
Last step. Type these two commands:
g.region rast=GRASS_DEM_FILEr.out.ascii input=GRASS_DEM_FILE output=DESIRED_NAME.grid null="-32747"
NOTES:
- GRASS_DEM_FILE - the name of the dem (as known by GRASS GIS) that you want to export.
- DESIRED_NAME - the name of the converted file - without the extension.
- -32747 - by default, any 'null' (unknown) values will be written as '*'. If you want somthing besides *, define it. In this example, we use -32747. Note, if * is fine with you, then don't add the null="XXX" part of the second command.
- The g.region command tells GRASS that the region of interest is the area contained by the GRASS_DEM_FILE. If you don't define this, the region will be assumed to be that of the MAPSET. Now, this may not be an issue, but, if the MAPSET covers a larger region than your 1/4 quad, you will have many, many more rows/columns of NULL values created by the r.out.ascii command. If the MAPSET and GRASS_DEM_FILE have the same boundaries, the command doesn't hurt.
SCENARIO 2 - Multiple LIDAR sets into one map set.
This section assumes you are familiar with BATCH PROCESSING With GRASS. If not, click on the previous link and take a quick look. Also, you should know how to create a MAPSET. If not, I'll be introducing it in the near future (hopefully).
1) Get the data files.
Essentially, follow the instructions seen in Scenario 1, Step 1. Only addition is a) you are downloading multiple LIDAR dems and b) you should keep track of which dems are on your north, south, east and west borders. If you don't do (b), then you are creating a lot more work for step 2.
As an aside, I'd probably create the following setup. Create a 'main' directory. For each Full quad, create a directory for that quad within the main directory. For each 1/4 dem, create a directory indicating is 1/4 quad (ne, nw, se, sw). Unzip the 1/4 dems in its respective 1/4 quad directory. This will help with the script example in a bit. So, for instance, for the NE 1/4 quad in the BreauxBridge Quad in St. Martin Parish, the directory would be StMartinParish/BreauxBridge/NE.
2) Gathering some information and importing into GRASS GISOnce again, we need a program call 'gdalinfo'. NOTE : commands will be italic and highlighted. NOTE: This is a command-line program.
Type: gdalinfo name_of_file.dem
You should have a lot of information casacading down the screen.
The items of importance are PIXEL_SIZE and Corner Coordinates. You will need those to create the GRASS location and MAPSET. If you don't know, please go HERE to find out how.What you are seeking to do is find the coordinates to create a rectangle region that incorporates ALL OF YOUR DEMS. Hence, you are looking for the northernmost, southermonst, eastmost and westmost points. This is why keeping track of the boundary DEMS (and which boundary the DEM is on) is useful; it cuts down the amount of work. Likely, a program could be written to do this; after all, GDAL is a library. However, I haven't gotten that far yet.
Once you have the coordinates, you can create the MAPSET. Note, I'd suggest adding some 'blank' space in the boundaries to the Mapset. At least for now. How much should be based on the PIXEL_SIZE. So, if your images are 5 by 5 meters, add 50 meters to the boundary. However, if your coordinate information is correct, you can just use the coordinates extracted.
Now, the fun part. You could type in the command used in Scenario 1, Step 2 several times. Or, you could use scripts. Two, to be exact. The first script actually generates the GRASS GIS script. The GRASS GIS script is then used in BATCH mode to import the data. (GENERATOR SCRIPT IS HERE)
So, assuming the following:
- we are going with the script approach.
- you created the suggested directory structure
- The script is in the directory containing the 'main' directory (so, in our example, the directory containing the 'StMartinParish' directory).
type the following command
sh constructInput.sh StMartinParish > input.sh
This will create the GRASS GIS script file. From here, you need to know
- The path to your Project and Location
- That you have read the BATCH MODE instructions
Run input.sh via batch mode. The data will now be imported!
3) Converting into friendly format.
Last step. Again, you could load up grass and run the commands in Scenario 1, Step 3 multiple time. Or you can run a couple of scripts. One to generate the GRASS GIS output script. The other being the GRASS GIS output script. NOTE: In this example, I'm assuming the GRASS Project name is stmartin and the MAPSET is PERMANENT
First, obtain the SCRIPT that generated the output script.
Second, run the following command
sh generateOutput.sh /path_to_grass_data_directory/stmartin/PERMANENT/cats > output.sh
where /path_to_grass_data_directory is the path to the GRASS GIS data (normally /something/grassdata).
Third, move the output.sh script to the directory where the ASCII files are to be placed.
Then, run output.sh in GRASS BATCH MODE.
And, you are done.
