I’ve never enjoyed working with contours. They seem to bog down my system more than any other layer type I work with. However, most of my clients are so used to looking at USGS Topo maps they expect to see them on at least one of the maps I produce for them. I recently worked on a project covering a five-town area in the Catskill Mountain region. The large area covered, and the ruggedness of the topography was proving exceptionally troublesome in processing their contours. So much so that I decided to look at other options to get the work done. I’ve used a variety of GIS tools over the years, but do most of my paying work exclusively in ArcGIS. It’s what I’m most familiar with, it does (nearly) everything I need it to do, and therefore provides my clients with the most efficient use of my time. However, in this situation that was not the case.
The one geoprocessing operation that frustrates me most often (in ArcGIS) is the Clip operation. It seems to take more time than most other geoprocessing tools, and often results in bad geometries. This happens so often, I usually resort to doing a union, and then deleting the unwanted areas of the Union results. For some reason this works much faster, and with more reliable results than doing a Clip.
Since what I wanted to do here was a clip on a contour layer, I was in for double trouble. Yes, I could have clipped the original DEM I wanted to produce the contours from first, then generated contours from the clipped DEM. But that wouldn’t have led to anything to write about. So, here’s a short comparison of how ArcGIS handled the process versus QGIS:
The hardware and software used:
ArcGIS 10, SP2
Windows 7, 64 bit
Dell Precision m2400 laptop
Intel Core 2 Duo CPU, 3.06GHz
8 GB RAM
QGIS 1.4.0
Ubuntu 11.4
Dell Inspiron 600m laptop
Intel Pentium M CPU, 1.60 Ghz
1GB RAM
A fair fight?
I started out with ArcGIS, and loaded up my 20’ contour lines and a 1 mile buffer of the study area to which I wanted to clip them. I began the clip operation 3 times. The first two times I had to cancel it because it was taking too long, and I needed to get some real work done. Curious to see how long it would really take, I let the process run overnight. The progress bar kept chugging away “Clip…Clip…Clip…Clip…”, and the Geoprocessing results window kept updating me with its progress, so I assumed it would complete eventually. In the morning, I looked in the Geoprocessing results window and found it had run for over 12 hours before throwing an error, never completing the clip operation. The error message said something about a bad geometry in the output. Really, no surprise there.
(Yes, those are lines in the picture above, not polygons. They’re very densely packed)
QGIS gets to play
The next day I decided to give QGIS a shot at it. I copied the two shapefiles over to my 6 year old lappy. (The contour.shp file was 1.3GB) fired up QGIS, and ran the Clip operation on the two files.
This time it took all of 17 minutes and 21 seconds to get a new contour layer.
So, who’s the winner here? Was it a far contest?
My take-away is, ESRI really needs to do some work on its Clip geoprocessing tool. As I said earlier, it is slow, and results in bad geometries more often than any of their other geoprocessing tools I use.
As ESRI is making it’s move away from VB Script and towards Python, I’m also slowly updating my stash of code snippets along the way. One of those little pieces of code I use quite often is one that identifies duplicate field names in a layer’s attribute table. I find this particularly helpful when I’m cleaning up tax parcel data, looking for duplicate parcel-ID numbers or SBL strings. Since I’ve been working a lot with parcel data lately, I figured it was time to move this code over to Python, too. So, here it is in step-by-step fashion…
1 – Add a new “Short Integer” type field to your attribute table (I usually call mine "Dup").
2 – Open the field calculator for the new field
3 – Choose "Python" as the Parser
4 – Check "Show Codeblock"
5 – Insert the following code into the "Pre-Logic Script Code:" text box making sure you preserve the indents:
uniqueList = [] def isDuplicate(inValue): if inValue in uniqueList: return 1 else: uniqueList.append(inValue) return 0
6 – In the lower text box, insert the following text, replacing "!InsertFieldToCheckHere!" with whichever field you want to check the duplicates for:
isDuplicate( !InsertFieldToCheckHere! )
This will populate your new field with 0′s and 1′s, with the 1′s identifying those fields that are duplicates.
One of the more popular analyses I’m asked to perform for my clients is a viewshed analysis. Beyond simply identifying what areas of a town are visible from roads or other public viewpoints, I’m often asked to help identify, and sometimes rank, areas that are most worthy of protection. One way to help a town identify and evaluate these high priority vistas, is to identify prominent ridgelines and the areas around them that are susceptible to inappropriate development.
One way to mitigate the impact of development on highly visible ridgelines is to make sure new buildings do not break the horizon line – the point where the ridge visibly meets the sky. Since most local zoning codes restrict building heights to 35-40 feet (in my client towns, anyway), producing a vertical buffer of 40 feet helps to identify the areas susceptible to such an intrusion.
The steps to produce such a vertical buffer are not overly complex, but I have not found them readily available online. So, for your benefit (and my easy reference) I outline the process here. (Note: I outline these steps specifically using ArcGIS 10. I’m sure it’s possible using other tools, but this is what I use most often in my daily work)
Begin with:
A DEM of your study area
This DEM must be an integer raster for one of the following steps, so start off by using the raster calculator.
(Arc Toolbox > Spatial Analyst Tools > Map Algebra > Raster Calculator)
Use the INT( “DEM” ) expression, where “DEM” is the elevation raster that you want to convert from a floating point to an integer raster.
Prominent ridgelines
Ridgelines are often generated using watershed boundaries, with extensive field checking to identify the more prominent features. For these calculations, the ridgelines must be in a raster format that includes the elevation of each raster cell. To convert a line-type ridgeline to a raster, first buffer the ridgeline to produce a polygon feature (I typically use a 20 foot buffer) and then use Extract by Mask.
(Arc Toolbox > Spatial Analyst Tools > Extraction > Extract by Mask)
Use the Integer DEM produced in the first step as the input raster, and the polygon buffer of the ridgelines for the feature mask data.
Generate the Euclidean Allocation of the Ridgeline Raster
This is where the need for an integer raster comes into play.
(Arc Toolbox > Spatial Analyst Tools > Distance > Euclidean Allocation)
Simply use the Ridgeline raster generated in the previous step as the input raster, and choose the elevation value as the Source field. This process will generate a new raster that covers the entire study area, with each cell holding the elevation value of the ridgeline raster that is nearest to it.
Generate the Vertical Buffer
Use the Raster Calculator to generate a vertical buffer of the ridgeline.
(Arc Toolbox > Spatial Analyst Tools > Map Algebra > Raster Calculator)
The expression will look something like this: “IntegerDEM” >= (“RidgelineEucoAllo”-12) where “Integer DEM” is the DEM produced in step 1, “RidgelineEucoAllo” is the Euclidean Allocation raster produced in step 3, and “12” is the height of the buffer you want to produce. In this case, the raster measurements are in meters, so using a buffer value of 12 results in about a 39 foot vertical buffer. This allows me to identify areas where there is the potential for a new building built to the maximum allowable height to break the horizon line.
Once the buffer is generated, there is usually some cleanup required. As you can see from these results, buffers are generated in areas not connected to the ridgeline areas we want to focus on, so I’ll delete these before moving on to the next steps in this town’s viewshed analysis.
I hope this is helpful, and as always, I welcome any comments and feedback.
Over the past week I’ve probably generated enough material for a half-dozen blog posts. However, since I have to get some billable hours in and invoices sent out this week, I’ll just post a short update for now on one aspect of project LOGAS. In my prior post I described the steps I took to get an Ubuntu based Apache/GeoServer up and running. Since then I’ve tweaked the process, and with a lot of help, been able to get everything in working order.
One of the roadblocks I faced was how to get Apache/Tomcat and GeoServer working together. Working individually was easy, together was not. As I progressed with my education in JavaScript, I ran into some cross-domain issues due to the fact my website and GeoServer were hosted on different machines. Even when I put them both on the same machine, the problems persisted due to the fact the webserver and the the map server use different ports (Apache – 80, and GeoServer – 8080). The solution turned out to be setting up a proxy server. Sounds easy. It was not.
After banging my head on the problem to the point of near concussion, a few online cohorts came to my rescue; in particular, @spara. She was generous enough to write up a script that installs Apache and the OpenGeo Suite, and configures Tomcat and the rest in a way that makes them all work together. She’s posted the code on github for anyone to use. Note: My previous post installed a full LAMP system. This one is limited to just Apache and the OpenGeo Suite, with some tweaks to the Tomcat servlet configuration.
I’ve tested her script on a few clean Ubuntu 10.10 setups, but I haven’t been able to get it to work as a single copy/paste yet. For me, it’s still a multi-step process. But, it does work and it is about as painless as it gets. What works consistently for me, is performing the OpenGeo Suite install first, and then running the Apache install/Tomcat configuration separately.
OpenGeo Suite
So far, the only process that works consistently for me is to first sudo to root, then setup the repository and install the OpenGeo suite. Start by typing the following into a terminal window:
After OpenGeo finishes installing, it will ask for a proxy URL (just leave it blank and hit enter unless you know what that means), a user name, and a password. Set these up as you wish. When all that’s done, move on the the next step.
Apache/Tomcat
Here’s the script for the Apache install/Tomcat setup, just copy and paste this into a terminal window:
sudo apt-get install -y apache2
sudo ln -s /etc/apache2/mods-available/proxy.conf /etc/apache2/mods-enabled/proxy.conf
sudo ln -s /etc/apache2/mods-available/proxy.load /etc/apache2/mods-enabled/proxy.load
sudo ln -s /etc/apache2/mods-available/proxy_http.load /etc/apache2/mods-enabled/proxy_http.load
sudo chmod 666 /etc/apache2/sites-available/default
sudo sed -i '$d' /etc/apache2/sites-available/default
sudo sh -c "echo ' ' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyRequests Off' >> /etc/apache2/sites-available/default"
sudo sh -c "echo '# Remember to turn the next line off if you are proxying to a NameVirtualHost' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPreserveHost On' >> /etc/apache2/sites-available/default"
sudo sh -c "echo ' ' >> /etc/apache2/sites-available/default"
sudo sh -c "echo '<Proxy *>' >> /etc/apache2/sites-available/default"
sudo sh -c "echo ' Order deny,allow' >> /etc/apache2/sites-available/default"
sudo sh -c "echo ' Allow from all' >> /etc/apache2/sites-available/default"
sudo sh -c "echo '</Proxy>' >> /etc/apache2/sites-available/default"
sudo sh -c "echo ' ' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /geoserver http://localhost:8080/geoserver' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /geoserver http://localhost:8080/geoserver' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /geoexplorer http://localhost:8080/geoexplorer' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /geoexplorer http://localhost:8080/geoexplorer' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /geoeditor http://localhost:8080/geoeditor' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /geoeditor http://localhost:8080/geoeditor' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /geowebcache http://localhost:8080/geowebcache' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /geowebcache http://localhost:8080/geowebcache' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /dashboard http://localhost:8080/dashboard' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /dashboard http://localhost:8080/dashboard' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPass /recipes http://localhost:8080/recipes' >> /etc/apache2/sites-available/default"
sudo sh -c "echo 'ProxyPassReverse /recipes http://localhost:8080/recipes' >> /etc/apache2/sites-available/default"
sudo sh -c "echo ' ' >> /etc/apache2/sites-available/default"
sudo sh -c "echo '</VirtualHost>' >> /etc/apache2/sites-available/default"
sudo chmod 644 /etc/apache2/sites-available/default
This installs Apache, and configures it and Tomcat so it appears as though your website and all of the OpenGeo Suite apps are using the same port. The last thing to do is test everything as I outlined in my previous post.
Test the Apache2 default website:
on the host computer – http://localhost – It works!
and via a remote computer (substitute your server’s domain here) – http://24.105.210.45/ – It works!
If everything is set up correctly, you should see something like this:
There are probably a few more tweaks that can be made to this process, but I don’t know enough about how Ubuntu works yet to make those changes. I do know that this process works most of the time. A restart between processes helps sometimes, and trying them in a different order occasionally works, too. If anyone can come up with a better, more streamlined approach, please feel free to let me know about it. For now, I’ll just have to live with a little complexity.
And in honor of the day, a little more cow bell, too -
I’ve run into a couple of roadblocks recently, regarding my experiments with my GeoSandbox. I want to be able to play with some of of the JavaScript libraries available, so I can continue my education on those fronts. Some cross domain issues arose when I started putting JavaScript in my web pages as my website is on a different server than is my GeoServer. So, that means learning some more about setting up web servers and how proxy servers work. Also, my old Dell 600m is starting to feel the effects of the increased demands placed on it as the GeoSandbox becomes more popular, and more complex. That means an eventual move onto a cloud server. Since cloud servers running an open source OS are much less expensive than those running Windows, I felt the need to begin transferring my entire setup over to open source tools.
For now, I’m just going to do a short outline of what it took to get a basic Ubuntu/Apache/OpenGeo Suite operating on my home server. Once I get a better handle on how Ubuntu and Apache work I’ll add some posts about that.
Starting with a clean install of Ubuntu 10.10
Download and install Ubuntu 10.10 Desktop
I downloaded Ubuntuu 10.10 Desktop Edition from their website:
Yes, I could have used the Server edition, but Coming from a Windows world, It’s easier for me to work with some kind of GUI than to go 100% command line. After downloading the .iso, I burned it to a CD, installed it on a clean hard drive, and then installed all updates.
Turn Ubuntu Desktop into an Apache Server
Install LAMP with a single command
In case you didn’t know, LAMP stands for “Linux Apache MySQL PHP”. This may or may not be more than I need for my purposes, but a one-line install looked like the easiest way to go, so:
Being somewhat new to Ubuntu, I missed what the first step under “Repository Setup” meant (sudo to root). Once I figured that out, things went smoothly. The steps are:
The last step was to see if everything was working properly, which it was.
Test the Apache2 default website:
on the host computer – http://localhost – It works!
and via a remote computer (substitute your server’s domain here) – http://24.105.210.45/ – It works!
In my case, everything worked as expected. After this, I continued trying to set up a custom website to serve my GeoServer maps, but ran into a few problems. I’ll be switching the GeoSandbox back and forth between the Windows and Apache server as I continue my climb up the learning curve, so don’t be surprised if some of the links appear to be broken on occasion.
