Mapping New York - georeferencing an old map

Last time, I mentioned that the Mohawk River State Park (former Schenectady Museum Nature Preserve) is missing from my map - because it appears neither in OpenStreetMap nor in the New York State inventory of public lands. How can I get it included?

Fortunately, I found online that someone had snapped a photograph of an old copy of the map that had once been posted at the Whitmyer Road parking lot, and posted it to a mountain bike forum. The picture looks to be output from a CAD system, or perhaps an early GIS. In any case, it appears to have been done to scale, and have just enough points that we can reference it to the map. Let's give it a try.



Pull up ‘Plugins->Georeferencer->Georeferencer’ from the QGIS menu bar.


Select the ‘Open Raster’ button in the toolbar, and load the ‘niskabike.jpg file that we downloaded from the mountain bike forum. Zoom in a couple of levels; you can use the ‘hand’ tool to pan around, just as you can with the map.

Now we have to identify corresponding points in the two maps. Click the ‘Add Point’ button in the toolbar, and click a point for which we can identify a corresponence. (I used the intersection of Whitmyer Drive and River Road as the first point.)

A dialog appears asking what corresponding point in the map to use. Click ‘From map canvas’ and select the point on the main map. The co-ordinates will be filled into the dialog box. Click ‘OK’

Continue this process for identifiable points that cover the map over a wide area. I used the turn of Whitmyer Road, the right angle in the upstream breakwater above the lock, the junction of the lock's retaining wall with the dam, the turn of Lock 7 Road near the parking area, and the junction of River Road and Rosendale Road.

Once you have a data table of corresponding points, you can let the georeferencer do its job by selecting ‘Start georeferencing’ in the toolbar.

You will need to fill out the dialog that appears, at least giving it a TIFF file for the output, and a projection to use (which should usually be the projection of the main project). The rest of the choices should be filled out as shown.

Clicking ‘OK’, we get the old map loaded up as a layer in the QGIS project. Lowering this layer in the stack (I put it below the contour lines), we can see that it aligns nearly perfectly. The retaining wall of the lock, the roads, and the bike path all superimpose nicely. The stream appears to have meandered a little, but it's pretty soggy around there, and I'll say that it's good enough.


Now we can start transcribing data by tracing over the old map. Create a new shapefile layer, choosing a Polygon layer, and adding a few attributes.

In the new layer, turn on edit mode (‘Toggle Editing’ in the toolbar). Create polygons that outline the land parcels of interest. Save the newly-created shapefile, and turn editing mode back off.

Lower the layer below the relief shading, style the land use appropriately, remove the temporary layer with the old map, and we're done.

(I could have transcribed the walking tracks, and so on, but I have other data sources for them. So read on to the next installment to see where that material is coming from.)

Read more...

Mapping New York - more polygons

Last time, we got the map to a pretty decent topographic road map. (In a future installment, I may discuss curating the OpenStreetMap data; it has recurrent errors that come from the Census Bureau TIGER files. But that's a side issue.) The next thing that I want to get into place is some polygon data: "who owns this land? What's it used for? Does the public have right of access? Where are the landmark buildings?

The relevant part of the first set of questions can mostly be answered by a database of publicly-owned lands. This database is getting into territory where the public databases aren't quite up to snuff. Within the Adirondack and Catskill Blue Line, the data are readily available from NYSGIS. The usual drill of using ‘ogr2ogr’ loads them into PostGIS:

ogr2ogr -f PostgreSQL -overwrite -t_srs EPSG:32618 \
    "PG:dbname=gis" DEC_Lands.shp \
    -nln nys_dec_lands -nlt MULTIPOLYGON -lco PRECISION=NO

Once again, the ‘-t_srs’ option is there to reproject the data into the projection that I intend to use for the finished map, and the ‘-lco PRECISION=NO’ works around a bug that causes a failure in inserting some of the numeric data.

Outside the Blue Line, the data come from a different place: the New York State Office of Cybersecurity. (I'd be intrigued to know why they became the custodian of the data.) In any case, they have a collection of files available on the NYSGIS web site.

Rather than putting each of these files into a separate table in the database, and hence needing a separate layer to show them, I decided to integrate them into a single table, and add a column to the data representing which data set a given row came from. For this, I decided to resort to scripting. Pulling out my handy-dandy Tcl interpreter, I ran the following loadall.tcl script:

set firsttime true
# Find all the shapefiles in the working directory
foreach file [glob *.shp] {

    # Extract the base name of each file, and create an 'ogr2ogr' command to load it
    set base [file rootname [file tail $file]]
    set cmd [list ogr2ogr -f "PostgreSQL"]
    if {$firsttime} {
        set firsttime false
        lappend cmd -overwrite
    } else {
        lappend cmd -append
    }
    # Add a 'data_source' column to identify which file we loaded
    lappend cmd -sql "SELECT *, '$base' AS data_source FROM $base" \
        -t_srs EPSG:32618 \
        -skipfailures \
        "PG:dbname=gis" \
        $file \
        -nln nys_public_land_boundaries \
        -nlt MULTIPOLYGON \
        -lco PRECISION=NO
    # Report on the console which file we're processing, and load it
    puts $cmd
    exec {*}$cmd >@stdout 2>@stderr
}

With these boundaries, what I mostly care about is “recreational” (go ahead and access) versus “nonrecreational” (permission needed, or special land use such as prisons and schools). (And I also want to treat the ‘AdirondackCatskill’ file specially, because that's the Blue Line, rather than reflecting public ownership.

Both of these layers need some styling. Rather than walk through that whole process, I have QML files attached at the end of the post.

Dealing with OpenStreetMap polygon data is rather more complicated, because it's got so many different things in the same file. I therefore made several different layers, with SQL queries to extract specific features.

(1) The first layer, I just left with the name, new_york_osm_polygon. This layer really represents "here are polygons that I don't know what to do with, yet." It is stacked behind everything else, and I usually leave it unchecked unless I'm actively working on styling for polygons. Its query looks like:

admin_level IS NULL
 AND ("boundary" IS NULL OR "boundary" NOT IN ('national_park'))
 AND ("waterway" IS NULL OR "waterway" IN ('boatyard','dam','dock','rapids','waterfall'))
 AND ("natural" IS NULL OR "natural" NOT IN ('bay', 'marsh', 'pond', 'swamp', 'water','waterway','wetland'))

This query excludes:

• Anything with an ‘admin_level’ attribute: these are administrative regions (states, counties, cities, towns, etc.)


• National park boundaries. I have these more accurately in the NYS Public Lands file.


• All waterways, other than man-made assets on the water.


• All wetlands, I've taken care of those already.

I style this letter either in bright yellow or icky purple, just to call attention to the unclassified features.

(2) The next layer is the layer where I describe land use. Its SQL query looks like:

"landuse" IS NOT NULL 
OR ("leisure" IS NOT NULL AND "leisure" NOT IN ('ice_rink','pitch', 'track', 'tennis_court'))
OR "aeroway" IS NOT NULL 
OR (amenity IN ('school', 'college', 'university', 'hospital') AND building IS NULL)

which translates to:

• Land use polygons


• Polygons marked 'leisure', except for a handful that usually appear inside parks and want to be rendered at a higher level.


• Polygons marked 'aeroway'.


• Polygons marked 'school', 'college', 'university', or 'hospital', except for buildings (these allow highlighting of campuses).

This set may need to be considered a work in progress; I expect these rules will need to be tweaked depending on the theme of the map.

For this layer, I created a fairly complex style with rules that fill the polygons in different colors according to land use. The QML is attached.

(3) Next up are the 'public lands' and 'DEC lands' layers, which I already discussed.

(4) Next, I have a few more types of region from OpenStreetMap. I call the layer OSM Subregion, and its query looks like:

(amenity IN ('parking')
OR leisure IN ('Dog Run', 'pitch' ,'tennis court','track'))
OR aeroway IN ('apron')
AND building IS NULL

so that it includes parking lots, dog runs, playing fields, tennis courts, racetracks, and airport aprons. What these types of object have in common is that they are usually layered atop another object (a shopping mall, industrial facility, park, airport, etc.), and so they look better rendered in an upper layer. Once again, I've attached a QML file to style them. This layer goes above the 'public lands' layers, but below the contour lines.

(5) Finally, there are polygons for a few man-made features:

building IS NOT NULL
OR leisure IN  ('pool','swimmin_pool','swimming_pool','wading_pool')

that is to say, buildings and swimming pools. This layer goes very high in the stack - above the roads. It provides footprints for these structures. Once again, I've attached the QML that styles it.

Wow, that's a lot of layers from one data set. But with them all in place, our map now has quite a lot of detail.

Next time, we'll make the map prettier, by adding shaded relief.

Attachments:

• NYS-dec-lands.qml - NYS DEC Lands


• NYS-public-land-boundaries.qml - NYS Public Land Boundaries


• OSM-building-style0.qml - OpenStreetMap - building footprints


• OSM-landuse-style0.qml - OpenStreetMap - land use polygons


• OSM-subregion-style0.qml - OpenStreetMap - subregions

Read more...