I came across an old backup set the other day and found a copy of the CVS -> Subversion repository I kept that included a lot of code that I wrote, inherited, or maintained.  The code is at least ten years old now so likely not of use to anyone.  I mainly did it to preserve the source for historical reasons.  If anyone is interested, you can find it at https://github.com/briangmaddox.

Yes, just in time for 2016… I’ve finished uploading and changed the links for the 2015 Census Tiger data converted to state or national datasets from the county files.  I’ve also stopped posting the single-file entries as there is not really a need to do that.  Hop over to https://brian.digitalmaddox.com/blog/?page_id=202 if you’re interested.

Mandy was in a car accident yesterday morning on her way home from work.  She’s ok, just got a bump on her head.  The car, however, didn’t fare so well, but hopefully with the insurance company it will be fixed and as good as new soon.

Hey all, I am finally in the process of updating my merged state- and national-level TIGER files to the 2014 data that they have put out.  You can find them at my GIS Data Page.  Note that Roads are not uploaded yet but I already updated the links on the download page so you will get 404 errors until I get them uploaded.  I cannot promise it will be tonight since I have to sleep sometime 😉  If you find any 404s on the others, let me know in case I missed a link.

As usual, these are my own value added files that I am publishing in case some people find them useful.  If you use these and your business fails, your wife leaves you, your dog dies, and you write a country music song about it, not my fault.

I’ll admit it, I really like old maps.  I especially like old topographic maps.  I think it started when I used to work for the US Geological Survey.  To me, it’s interesting to see how things change over time.  From my old urban grown prediction days, I like to see where and when populations change.

Since the USGS put out their Historical Topographic Map Collection (HTMC), I’ve been playing with the maps off and on for about a year now.  I finally decided to start merging available maps of the same scale and vintage to study and possibly do feature extraction down the road.  I’ll be placing them for download here in case anyone is interested as I process them.

I thought I’d share how I put them together in case anyone else is interested.  The first step is to go to the USGS website and pick the files you want to use.  The files there are available in GeoPDF format.  First thing you need to understand is that you may not find a map covering your area of interest at your scale of interest and vintage.  Not everything managed to survive to the current day.  For example, I made a merged 125K map of Virginia and most of southern VA is missing at that resolution.

Once I download the GeoPDFs I want to work on, I use a modified version of the geopdf2gtiff.pl script from the Xastir project.  The link to my modifications can be found here.  I use LZW compression for my GeoTIFFs as it’s lossless and keeps the quality from the GeoPDFs.  It is a Perl script and requires that you have GDAL and the GDAL API installed.  What it does is calculate the neat-lines of the GeoPDF and then clips it to the neat-line while converting it to a GeoTIFF.  Running it as as simple as:

geopdf2gtiff.pl inputfile.pdf

Once you have all of your GeoPDF files download and converted, the next step is to merge them.  The fastest way I’ve found to merge involves using gdalbuildvrt and gdal_translate, also from GDAL.  The first step is to create a virtual dataset of all of your files by running something like:

gdalbuildvrt -resolution highest -a_srs EPSG:4326 merged.vrt parts/*.tif

The options I chose here are to pick the highest pixel resolution (-resolution) based on the input files.  For this case the resolutions should be the same, but this way I don’t have to go through and verify that.  Next I change the projection of the output file to WGS84 (-a_srs).  Next is the file name of the virtual dataset and then the input files.

Now that the virtual dataset is done, it’s time to actually merge all of the files together.  The virtual dataset contains the calculated bounding box that will contain all of the input files.  Now we use gdal_translate to actually create the merged GeoTIFF file:

gdal_translate -of GTiff -co COMPRESS=LZW -co PREDICTOR=2 merged.vrt ~/merged.tif

Here again I use LZW compression to losslessly compress the output data.  Note that gdal_translate will automatically add an Alpha channel as Band 4 in the image to denote areas that had no input data.  That’s why we do NOT add the -addalpha flag to gdalbuildvrt.  For performance tips, I’d suggest keeping the source data and output file on separate drives unless you’re running something like a solid state drive.  To give you an idea of the output file sizes, Virginia merged (which did have a lot of areas missing), was around 500 megabytes.

Next you’ll need a Shapefile to use as a cut file to clip the data.  Since I have the Census Tiger 2013 data in a local PostGIS database (see previous posts to this blog), I used QGIS to select just the VA state outline and then saved it as a Shapefile.

Finally, we will use gdalwarp to clip the merged GeoTIFF against the state outline to produce the clipped GeoTIFF that is just the state itself.  This operation can take a bit of time depending on how powerful a machine you’re running it on.  The command you will use is similar to this:

gdalwarp --config GDAL_CACHEMAX 1024 -wm 1024 -cutline va_outline.shp -crop_to_cutline -multi -t_srs EPSG:4326 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=YES -co TILED=YES ~/merged.tif clipped.tif

Some of the command line parameters I used are optional, I just tend to leave them in since I do a lot of copying and pasting 😉  First we tell GDAL to increase the size of its caches using the –config GDAL_CACHEMAX and -wm options.  Next we specify the file to clip against with the -cutline and -crop_to_cutline options.  The -multi option tells GDAL to process using multiple threads.  I again specify the output projection and the LZW compression parameters.  Here I also specify the BIGTIFF option just in case the output file goes over four gigabytes.  Finally, I tell gdalwarp to tile the output TIFF so it will load faster in a GIS by separating it into tiles.

The output will look something like the below figure.  I’ll start posting files as I get time.  Hope everyone is having a great holiday!

Clipped Virginia 125K Old Maps

I’ve been keeping up with tracking how many botnets are out there scanning WordPress blogs.  I’ve eventually resorted to blocking huge chunks of the Internet via .htaccess files.   So far it’s been quite effective in limiting the number of botnet login attempts.

If anyone is interested, I’ve put the limit portion of my .htaccess file here.  Feel free to use it and customize for your needs.

Lately I’ve been getting tons and tons of login attempts from what appear to be botnets.  Since I’m getting tired of banning the IPs individually, I’m temporary taking to banning entire countries and ISP’s from hitting my blog.  If you’re in that group, sorry guys.  Take it up with your ISP.

Here are some stats I’ve been gathering.

IP Addresses Grouped by ISPs

As you can see, I get a bunch from Godaddy and French ISP Ovh.  I’ve also banned Godaddy IP’s, Ovh, and Media Temple.  I’ll be adding others once I find all of their allocated net ranges.

For reference, here’s a copy of my current list along with attempts:

Guess I should be flattered that I’m getting all of this “attention” 🙂

Towards the end of my tenure at the US Geological Survey, I was the project manager and principal investigator of Restoration of Data from Lossy Compression.  The goal of the project was to find ways to restore fine detail that was lost during lossy compression processes such as JPEG.  I had submitted an Open File report through the review process, but left the USGS in 2006 before the paper had completed review.  As I had left, it basically fell through the cracks and was never officially published.

I had forgotten about it until recently when I was updating my resume.  So, without further ado, I have put the paper here.  I took out the USGS logo and what not since it was never officially published by them.  So for a flashback into what I was doing in 2006, have fun reading it 🙂

Before going into what tools and data are available, I feel a quick introduction to geospatial data is in order. I’m going to gloss over some things and greatly oversimplify others. If you want to know more, I would suggest checking out books such as Map Projections: A Working Manual by Snyder or some excellent tutorials on map projections and GIS that can be found via your favorite web search engine.

To begin, there are two main types of geospatial data available: raster and vector. Raster data are usually things such as aerial or satellite photographs, scanned paper maps, and so on. Vector data are mathematical vectors that describe items such as lines for roads, polygons for areas like parks, and others. As is, aerial photographs or vectors do not fully describe the Earth’s surface. They do not represent a one to one correlation between a point on the map and a location on the ground. They have to be run through a process called georectification to mathematically translate them to a coordinate system known as a map projection.

USGS DOQQ of the Washington DC Area Courtesy of Wikipedia

Map projections are mathematical descriptions of the surface of the Earth. As most of you (hopefully) know, the Earth is not flat. It’s a large spherical object with dents (canyons) and bumps (mountains). Since it’s not a flat sphere, it’s impossible to model it with 100 percent accuracy. Map projections were created that describe the Earth in terms of a series of equations that define a spheroid with various characteristics. These characteristics are useful for different types of representations, from maps of continents to maps of cities. I suggest you look at websites such as Map Projections: From Spherical Earth to Flat Map to get a more in-depth discussion. Basically, think of a map projection as establishing an average elevation of the Earth’s surface with areas where it’s more accurate than others.

USGS Traverse Mercator Projection. Courtesy Wikipedia

Georectification then takes the flat two-dimensional image or vector data and warps it so that it fits the three-dimensional surface of the Earth as approximated by the map projection. Once transformed, each pixel in a map image or each point on a vector matches its corresponding point on the Earth’s surface. Inside a GIS, you can then move the cursor around and see the matching latitude and longitude on the ground. If you wanted, you could take a GPS and even go there.

Axonometric Projection of a 3D object to a 2D surface.  Courtesy Wikipedia

GIS file formats then allow the mathematical parameters that went into georectification to be encoded into the file as well. These are usually a large amount of numbers that describe the equations used and the inputs into all of those equations. Many times you’ll see something like below if you use a tool such as listgeo that can parse the information out of the file format. Note that most formats assume a Cartesian coordinate system running from coordinates (0,0) in the upper left corner to some X and Y to the size of the data set.

Cartesian Coordinate System Courtesy Wikipedia

listgeo o39102g6.tif
Geotiff_Information:
Version: 1
Key_Revision: 0.2
Tagged_Information:
ModelTiepointTag (2,3):
0 0 0 
691331.977 4417194.85 0 
ModelPixelScaleTag (1,3):
2.4384 2.4384 0 
End_Of_Tags.
Keyed_Information:
GTModelTypeGeoKey (Short,1): ModelTypeProjected
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
ProjectedCSTypeGeoKey (Short,1): PCS_NAD27_UTM_zone_13N
PCSCitationGeoKey (Ascii,25): "UTM Zone 13 N with NAD27"
End_Of_Keys.
End_Of_Geotiff.

PCS = 26713 (NAD27 / UTM zone 13N)
Projection = 16013 (UTM zone 13N)
Projection Method: CT_TransverseMercator
ProjNatOriginLatGeoKey: 0.000000 ( 0d 0' 0.00"N)
ProjNatOriginLongGeoKey: -105.000000 (105d 0' 0.00"W)
ProjScaleAtNatOriginGeoKey: 0.999600
ProjFalseEastingGeoKey: 500000.000000 m
ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4267/NAD27
Datum: 6267/North American Datum 1927
Ellipsoid: 7008/Clarke 1866 (6378206.40,6356583.80)
Prime Meridian: 8901/Greenwich (0.000000/ 0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left ( 691331.977, 4417194.851) (102d45'44.74"W, 39d53' 6.48"N)
Lower Left ( 691331.977, 4400540.579) (102d46' 2.23"W, 39d44' 6.68"N)
Upper Right ( 704806.576, 4417194.851) (102d36'17.89"W, 39d52'55.15"N)
Lower Right ( 704806.576, 4400540.579) (102d36'36.61"W, 39d43'55.41"N)
Center ( 698069.276, 4408867.715) (102d41'10.37"W, 39d48'31.03"N)

From the above output, we can see a few things that describe this image. It came from a USGS DRG that was mathematically warped to a projection (which is why some of them look “tilted” if you’ve looked at a lot of DRGs). Under the ModelTiepointTag key we see several zeros then the numbers 691331.977 and 4417194.85. If you look at the DRG as a two dimensional grid, those numbers represent the latitude and longitude (x, y) coordinates of the upper left pixel in the image. The projection of the DRGs uses meters so the x and y here represent meters in the coordinate system. Under the ModelPixelScaleTag we see another set of numbers, including 2.4384 listed twice. This is the per-pixel scale of the image in meters. So this means that as you increase in the x- or y- direction of the image, each pixel you count up equals 2.4384 meters on the ground. Moving from pixel (0,0) to (10,10) in this case only moves you by a value of 10 on the grid, but would move you 24.384 meters on the ground.

Scale is another issue that may seem counter-intuitive to many people. There are two scales of data that you may hear about: large and small scale. Large scale actually means data that is more detailed than a small scale. The scales do not reference the size of the map but the scale of the map to the Earth. Obviously a 1:1 scale map of the Earth is not feasible since it would have to be as large as the planet. Maps are scaled down to make them much smaller. A common size of a large scale map is 1:25,000 and a small scale map is 1:100,000. The 1:25,000 map is called a large scale because it’s fractional size is larger than that of the 1:100,000 scale map (think of the scales as fractions). A 1:25,000 large scale map does not show as much surface area as the 1:100,000 small scale map, but it shows an area in much more detail than the small scale map does. Think of a large scale map as maps of things like cities, while small scale maps are more along the lines of maps of entire countries.

With the scale discussion out-of-the-way, we then must turn to the usability of geospatial data. Much digital map data was collected from paper maps. Thus, the digital data has a scale where it is most accurate and less accurate at others. Data collected at 1:25,000 scale will be more accurate than data collected at 1:100,000. A GIS allows you to zoom in to the data and in some cases over-zoom past the actual accuracy of the data. Consider data such as a USGS Digital Elevation Model that was collected on a ten meter grid. You can zoom to a view where each pixel on the screen is smaller than ten meters on the ground, but your measurements will not be any more accurate than the ten meter grid you are viewing. Another example is looking at two road maps where one was collected at 1:25,000 and the other at 1:250,000. If you overlay the road maps on to a georeferenced aerial photograph, the 1:25,000 lines will more closely follow actual roads than the ones from the 1:250,000. You cannot expect the 1:250,000 data to be as accurate as the larger scale data, so zooming into a data set will not make the 1:250,000 data any more accurate than the scale where it was captured.

With this out-of-the-way, next time we will start looking at Open Source tools to use the data.