1<!DOCTYPE html> 2<html lang="en"> 3 <head> 4 <meta property="fb:pages" content="520981908088655" /> 5 <meta charset="utf-8"> 6 <title>My Musings</title> 7 <meta name="viewport" content="width=device-width, initial-scale=1.0"> 8 <meta name="description" content=""> 9 <meta name="author" content=""> 10 11 <!-- Le styles --> 12 <link href="/css/bootstrap-spacelab.min.css" rel="stylesheet"> 13 <style> 14 body { 15 padding-top: 60px; /* 60px to make the container go all the way to the bottom of the topbar */ 16 } 17 </style> 18 <link href="/css/bootstrap-responsive.css" rel="stylesheet"> 19 <link href="/css/custom.css" rel="stylesheet"> 20 21 <!-- HTML5 shim, for IE6-8 support of HTML5 elements --> 22 <!--[if lt IE 9]> 23
23<script src="http://html5shim.googlecode.com/svn/trunk/html5.js"></script>
23 24 <![endif]--> 25 26 <!-- Fav and touch icons --> 27 <link rel="apple-touch-icon-precomposed" sizes="144x144" href="/ico/apple-touch-icon-144-precomposed.png"> 28 <link rel="apple-touch-icon-precomposed" sizes="114x114" href="/ico/apple-touch-icon-114-precomposed.png"> 29 <link rel="apple-touch-icon-precomposed" sizes="72x72" href="/ico/apple-touch-icon-72-precomposed.png"> 30 <link rel="apple-touch-icon-precomposed" href="/ico/apple-touch-icon-57-precomposed.png"> 31 <link rel="shortcut icon" href="/ico/favicon.png"> 32 <link rel="stylesheet" href="/css/prism.css"> 33
33<script type="text/javascript" async src="/js/prism.js"></script>
vendor: 1 bytes, line 33
33
34<script type="text/javascript" async 35 src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-MML-AM_CHTML"> 36</script>
36 37
37<script src="/js/code.jquery.com/jquery-1.12.3.min.js"></script>
37 38
38<script src="/js/cdn.datatables.net/1.10.12/js/jquery.dataTables.min.js"></script>
38 39 <link rel="stylesheet" href="/css/cdn.datatables.net/1.10.12/css/jquery.dataTables.min.css" /> 40
40<script> 41 $(function(){ 42 $(".datatable").DataTable({"pageLength": 50}); 43 $(".layerck").change(function(e) { 44 layers_dict[$(this).attr('id')].setVisible($(this).is(":checked")); 45 }); 46 }); 47 </script>
47 48 49 50 51 </head> 52 53 <body> 54 <!-- Google Tag Manager -->
vendor: 66 bytes, lines 54-56
54 55 <noscript><iframe 56 src="//www.googletagmanager.com/ns.html?id=
56GTM-KKPT75
vendor: 89 bytes, lines 56-59
56" 57 height="0" width="0" 58 style="display:none;visibility:hidden"></iframe></noscript> 59
59<script>
vendor: 331 bytes, lines 59-64
59(function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start': 60 new Date().getTime(),event:'gtm.js'});var 61 f=d.getElementsByTagName(s)[0], 62 j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src= 63 '//www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f); 64 })(window,document,'script','dataLayer','
64GTM-KKPT75
vendor: 3 bytes, line 64
64');
64</script>
64 65 <!-- End Google Tag Manager --> 66 67 <div class="navbar navbar-inverse navbar-fixed-top"> 68 <div class="navbar-inner"> 69 <div class="container"> 70 <a class="btn btn-navbar" data-toggle="collapse" data-target=".nav-collapse"> 71 <span class="icon-bar"></span> 72 <span class="icon-bar"></span> 73 <span class="icon-bar"></span> 74 </a> 75 <a class="brand" href="/">My Musings</a> 76 <div class="nav-collapse collapse"> 77 <ul class="nav"> 78 <li class="active"><a href="/">Home</a></li> 79 <li class=""><a href="/posts">Posts</a></li> 80 <li class=""><a href="/about">About</a></li> 81 <li class=""><a href="/talks-given">Talks</a></li> 82 </ul> 83 </div><!--/.nav-collapse --> 84 </div> 85 </div> 86 </div> 87 88 <div class="container"> 89 90 <div class="row"> 91 <div class="span2"> 92 <h3>About Me!</h1> 93 <img src="/img/me2.jpg"> 94 <p class="mini-about-me">This blog is about my musings and thoughts. I hope you find it useful, at most, and entertaining, at least.</p> 95 <p class="resume"><a href="/resume">Résumé</a> [<a href="/pdf/James_Keener_Resume.pdf">PDF</a>]</a></p> 96 <h4>Other Pages</h4> 97 <p><a href="/quotes">Quotes</a></p> 98 <p><a href="/links">Links</a></p> 99 <h4>Presence Elsewhere</h4> 100 <p><a href="mailto:jim@jimkeener.com">jim@jimkeener.com</a></p> 101 <p><a href="https://www.github.com/jimktrains">GitHub</a></p> 102 <p><a href="https://bitbucket.org/jimktrains">BitBucket</a></p> 103 </div> 104 <div class="span8"> 105 <h1><a href="/posts/luddite2">Why I'm a Luddite, kind-of (Part 2)</a></h1> 106<p>Date: <a href="/posts/by-date/2022">2022</a><del><a href="/posts/by-date/2022/05">05</a></del><a href="/posts/by-date/2022/05/27">27</a><br /> 107<br /> 108Tags: <a href="/posts/tags/technology">technology</a> <a href="/posts/tags/Luddite">Luddite</a> <a href="/posts/tags/freedom">freedom</a> <a href="/posts/tags/liberty">liberty</a></p> 109 110<p>In 2013 I wrote a ~~ramble~~ post called <a href="/posts/luddite">Why I’m a Luddite, 111kind-of</a>
111. Now, some 9 years later, I want to revisit<br /> 112some of my opinions then and add some more color to my basic premise. If<br /> 113anything, I think my opinion has gone from “sometimes technology doesn’t<br /> 114actually make things better” to “usually technology doesn’t actually<br /> 115make things better.”</p><br /> 116<p>I’ve come to loathe “smart” interfaces.</p><br /> 117<p>In 2013 I was still two years away from having a child. While Anne and I<br /> 118weren’t able to get our families to strictly adhere to the “no batteries<br /> 119allowed” rule, they did follow it pretty well and I’m pretty happy with<br /> 120the result. Our children watch a little more TV than I’d prefer. A<br /> 121significant portion of it is made up of history and science<br /> 122documentaries and they are (usually) more than willing or even beg to go<br /> 123outside to play, so it’s a mixed bag.</p><br /> 124<p>One sentence, of many, I do think has aged very poorly and doesn’t really reflect<br /> 125my current beliefs, though, is “I just wish so many people wouldn’t use<br /> 126them as a substitute for real life interaction.” Before 2020, but<br /> 127specially in a post-2020 World, that sentence was just me talking out of<br /> 128my butt. While I still think real-life interactions are important,<br /> 129interacting with others any way in which we’re comfortable is a<br /> 130positive.</p> 131 132<h1><a href="/posts/visualizing-system-wide-headway">Visualizing System-wide Headway</a></h1> 133<p>Date: <a href="/posts/by-date/2019">2019</a><del><a href="/posts/by-date/2019/11">11</a></del><a href="/posts/by-date/2019/11/13">13</a><br /> 134<br /> 135Tags: <a href="/posts/tags/postgis">postgis</a> <a href="/posts/tags/postgres">postgres</a> <a href="/posts/tags/gis">gis</a> <a href="/posts/tags/gtfs">gtfs</a> <a href="/posts/tags/transit">transit</a></p> 136 137 <p><p>I was bored this evening and decided to visualized the system-wide headway,<br /> 138that is the amount of time between buses at a stop for every stop in the<br /> 139system.</p><br /> 140<p>Given the <a href="/posts/stop-spacing-part-1">stop-distance table</a> from the previous<br /> 141post, I set out to compute the headway for each stop segment (the portion<br /> 142of a route between stops).</p><br /> 143<pre><code>— Compute the max headway for each stop per hour.<br /> 144create table overall_headway as <br /> 145select stop_id, 146 hour, 147 max(headway) as max_headway<br /> 148from ( 149 select stop_id, 150 extract(hour from departure_time) as hour, 151 departure_time, 152 lag(departure_time, 1) over w as prev_departure_time, 153 extract(epoch from (departure_time – lag(departure_time, 1) over w))/60 as headway 154 from gtfs.stops_3365 155 join gtfs.stop_times using (stop_id) 156 window w as (partition by stop_id order by departure_time asc)<br /> 157)x <br /> 158group by stop_id, hour;</p> 159 160 161 <p>— Assign the segments created for the stop-distance calculation<br /> 162— a headway from overall_headway.<br /> 163create table stop_line_headway as<br /> 164select coalesce(ohs.hour, ohp.hour) as hour, 165 sdl.stop_id, 166 sdl.dist, 167 prev_stop_id, 168 min(least(ohs.max_headway, ohp.max_headway)) as headway, 169 the_geom<br /> 170from stop_dist_line sdl<br /> 171left join overall_headway ohs 172 on sdl.stop_id = ohs.stop_id<br /> 173left join overall_headway ohp 174 on sdl.prev_stop_id = ohp.stop_id<br /> 175where (ohs.hour is null or ohp.hour is null or (ohs.hour = ohp.hour)) 176 and st_length(the_geom) < 2500<br /> 177group by coalesce(ohs.hour, ohp.hour), sdl.stop_id, sdl.dist, prev_stop_id, the_geom;</p> 178 179 <p>alter table stop_line_headway add id serial primary key;<br /> 180</code></pre><br /> 181<p>I then pulled this into <span class="caps">QGIS</span>, colored by graduation, and filtered per hour.<br /> 182For each our, I exported a <span class="caps">PNG</span> of the map and tried by hand at <span class="caps">KDE</span>nlive for the<br /> 183first time. (Although <code>ffmpeg</code> could have gotten the job done as well.)</p><br /> 184<iframe width="560" height="315" src="https://www.youtube.com/embed/hZrRQW25Z_0" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe></p> 185 186<p>Raw video, images, and more details on building the map to come.</p> 187 188<h1><a href="/posts/stop-spacing-part-1">Stop Spacing, Part 1</a></h1> 189<p>Date: <a href="/posts/by-date/2019">2019</a><del><a href="/posts/by-date/2019/11">11</a></del><a href="/posts/by-date/2019/11/12">12</a><br /> 190<br /> 191Tags: <a href="/posts/tags/postgis">postgis</a> <a href="/posts/tags/postgres">postgres</a> <a href="/posts/tags/gis">gis</a> <a href="/posts/tags/gtfs">gtfs</a> <a href="/posts/tags/transit">transit</a></p> 192 193 <p><p>For reasons I’ll get into in another post, I wanted to have a rough idea of the<br /> 194distribution of the distance between stops. To do this I used the<br /> 195<a href="https://gtfs.org/"><span class="caps">GTFS</span></a> import for the <a href="https://www.portauthority.org/business-center/developer-resources/">Port Authority of Allegheny County<br /> 196GTFS</a> that<br /> 197I did earlier this year when <a href="/posts/calculating-service-area-part-1-introduction">calculating the service<br /> 198area</a>.</p><br /> 199<p>A quick note about projections. Since I had created a table with a geometry in<br /> 200the Pennsylvania South (ftUS) state plane( <a href="https://epsg.io/3365"><span class="caps">EPGS</span> 3365</a>),<br /> 201any distance calculations I do will be done in feet.</p><br /> 202<p>A quick, rough and dirty experiment went like this:</p><br /> 203<pre><code>create table stop_dist_line as <br /> 204select stop.stop_id, 205 stop.stop_name, 206 prev.stop_id as prev_stop_id, 207 prev.stop_name as prev_stop_name, 208 — Just using euclidean distance right now. 209 — This could be improved by using the on-the-street distance. 210 stop.the_geom <-> prev.the_geom as dist, 211 st_setsrid(st_makeline(stop.the_geom, prev.the_geom), 3365) as the_geom<br /> 212from ( 213 — This gets unique pairs of stops, regardless of ordering. 214 select distinct least(stop_id, prev_stop_id) as prev_stop_id, 215 greatest(stop_id, prev_stop_id) as stop_id 216 — This gets all pairs of stops. In theory there could be duplicate pairs 217 — based on direction. 218 from (select stop_id, 219 lag(stop_id, 1) over stop_order prev_stop_id 220 from gtfs.stop_times 221 join gtfs.trips using (trip_id) 222 window stop_order as (partition by trip_id, direction_id order by stop_sequence) 223 ) unordered_pairs<br /> 224) pairs<br /> 225join gtfs.stops_3365 stop on stop.stop_id = pairs.stop_id<br /> 226join gtfs.stops_3365 prev on prev.stop_id = pairs.prev_stop_id;</p> 227 228 <p>— Remove segments that don’t have a line because there is no<br /> 229— previous stop to have a line or have a distance from.<br /> 230delete from stop_dist_line where prev_stop_id is null;</p> 231 232 <p>— Register the new geometry column<br /> 233select Populate_Geometry_Columns();<br /> 234— Add a PK to play more nicely with QGIS<br /> 235alter table stop_dist_line add id serial primary key;<br /> 236— Add an index on the geometry.<br /> 237create index on stop_dist_line using gist(the_geom);<br /> 238</code></pre><br /> 239<p>We can then get some idea of the distribution by bucketing by 100s.</p><br /> 240<pre><code>select dist_group, 241 ttlsf, 242 ttl, 243 round(ttlsf/ttl, 2) as perc, 244 c 245from 246 (select dist_group, 247 sum(c) over (rows between unbounded preceding and current row) as ttlsf, 248 sum(c) over () as ttl, 249 c 250 from ( 251 select 100*floor(dist / 100) as dist_group, 252 count(*) as c 253 from stop_dist_line 254 group by dist_group 255 order by dist_group) x)y; 256</code></pre><br /> 257<p>Giving</p><br /> 258<pre><code> dist_group | ttlsf | ttl | perc | c 259------------+-------+------+------+------ 260 0 | 149 | 7732 | 0.02 | 149 261 100 | 223 | 7732 | 0.03 | 74
262 200 | 681 | 7732 | 0.09 | 458 263 300 | 1568 | 7732 | 0.20 | 887 264 400 | 2528 | 7732 | 0.33 | 960 265 500 | 3579 | 7732 | 0.46 | 1051 266 600 | 4462 | 7732 | 0.58 | 883 267 700 | 5135 | 7732 | 0.66 | 673 268 800 | 5589 | 7732 | 0.72 | 454 269 900 | 5946 | 7732 | 0.77 | 357 270 1000 | 6234 | 7732 | 0.81 | 288 271 1100 | 6456 | 7732 | 0.83 | 222 272 1200 | 6628 | 7732 | 0.86 | 172 273 1300 | 6766 | 7732 | 0.88 | 138 274 1400 | 6852 | 7732 | 0.89 | 86 275 1500 | 6921 | 7732 | 0.90 | 69 276 ... 277</code></pre><br /> 278<p>However, because PostgreSQL is awesome, we can also just compute percentiles<br /> 279directly.</p><br /> 280<pre><code>select percentile_cont(0.25) within group (order by dist) as p_25, 281 percentile_cont(0.35) within group (order by dist) as p_35, 282 percentile_cont(0.45) within group (order by dist) as p_45, 283 percentile_cont(0.50) within group (order by dist) as p_50 284from stop_dist_line; 285</code></pre><br /> 286<p>Giving</p><br /> 287<pre><code> p_25 | p_35 | p_45 | p_50 288------------------+------------------+------------------+------------------ 289 438.775035648918 | 518.085656645151 | 590.352075165022 | 631.007847560984 290</code></pre><br /> 291<p>Next steps would be to correct many of the little corners cut, e.g. not using<br /> 292on-the-street distance, and rerunning the calculations.</p></p> 293 294<h1><a href="/posts/calculating-service-area-part-4-next-steps">Calculating Service Area, Part 4: Next Steps</a></h1> 295<p>Date: <a href="/posts/by-date/2019">2019</a><del><a href="/posts/by-date/2019/01">01</a></del><a href="/posts/by-date/2019/01/27">27</a><br /> 296<br /> 297Tags: <a href="/posts/tags/postgis">postgis</a> <a href="/posts/tags/postgres">postgres</a> <a href="/posts/tags/transit">transit</a> <a href="/posts/tags/census">census</a></p> 298 299<p>In <a href="/posts/calculating-service-area-part-1-introduction">Part 1</a> we<br /> 300talked about the project and got some of the tools and data we need to<br /> 301start. A quick reminder of our motivation: to figure out the percentage<br /> 302area and population the county transit agency serves by using a<br /> 303“walkshed”. </p><br /> 304<h2>Parts</h2><br /> 305<ul><br /> 306<li><a href="/posts/calculating-service-area-part-1-introduction">Part 1: Introduction and Setup</a></li><br /> 307<li><a href="/posts/calculating-service-area-part-2-defining-the-walkshed">Part 2: Defining the Walkshed</a></li><br /> 308<li><a href="/posts/calculating-service-area-part-3-demographics">Part 3: Demographics</a></li><br /> 309<li>Part 4: Next Steps</li><br /> 310</ul><br /> 311<h2>Walkshed</h2><br /> 312<p>There are two changes I should make to the walkshed calculation before I dig<br /> 313into per-route statistics: (1) excluding interstates and (2) light rail and<br /> 314busways stations should have a half mile walkshed. I don’t believe they made a<br /> 315difference at the county level; once I redo the calculations I’ll come back to<br /> 316see if I have to eat those words.</p><br /> 317<h2>More In-Depth Demographic Analysis</h2><br /> 318<p>Since we have all of the census data, all of the route data, and we broke our<br /> 319walkshed down by stop it’d be quite easy to look at the service area as a<br /> 320whole, and routes in particular against other socio-economic indicators.</p><br /> 321<p>Also, more interestingly, who <em>isn’t</em> being served?</p><br /> 322<h2>Filter Walkshed by Reliable Transit</h2><br /> 323<p>I’d love to look at what the walkshed is when only stops that are served more<br /> 324frequently than every 30 min mid-day is included. Another alternative would be<br /> 325late at night, or weekends. Essentially, what’s the service area that you could<br /> 326do a good deal in without having to get into a car?</p><br /> 327<h2>Build a web-based tool to export this data</h2><br /> 328<p>I would love to make this data more usable and accessible. I previously made a<br /> 329<a href="https://github.com/jimktrains/pptmap">little tool</a> using python, mapnik, and a<br /> 330file from GDAL2Tiles. I’d like to explore if <span class="caps">QGIS</span>’ “save project to database“<br /> 331functionality could be used (<em>i.e.</em> mangled, forced, hacked?) to do this more<br /> 332cleanly.</p><br /> 333<h2>Isochron / Isotrip maps</h2><br /> 334<p>I’ve been wanting to build an application where you could pick a location<br /> 335and it’d generate two (interactive) maps of all other stops: an isochron map<br /> 336showing how long it’d take to get anywhere else (and that has 15-min isochron<br /> 337lines showing the boundry of how long it takes to get how far) and an isotrip<br /> 338map showing how many trips/transfers it’d take to get to the other stop, with<br />
339similar isotrip lines.</p><br /> 340<h2>Create Subject Tables</h2><br /> 341<p>The Access DBs also contain information on what all the fields mean and which<br /> 342subject tables they belong to. We can extract it with the <code>mdb-export</code> tool.</p> 343<pre><code>mdb-export SF1_Access2003.mdb data_field_descriptors &gt; data_field_descriptors.csv 344</code></pre> 345 346<h1><a href="/posts/calculating-service-area-part-3-demographics">Calculating Service Area, Part 3: Demographics</a></h1> 347<p>Date: <a href="/posts/by-date/2019">2019</a><del><a href="/posts/by-date/2019/01">01</a></del><a href="/posts/by-date/2019/01/27">27</a><br /> 348<br /> 349Tags: <a href="/posts/tags/postgis">postgis</a> <a href="/posts/tags/postgres">postgres</a> <a href="/posts/tags/transit">transit</a> <a href="/posts/tags/census">census</a></p> 350 351 <p><p>In <a href="/posts/calculating-service-area-part-1-introduction">Part 1</a> we<br /> 352talked about the project and got some of the tools and data we need to<br /> 353start. A quick reminder of our motivation: to figure out the percentage<br /> 354area and population the county transit agency serves by using a<br /> 355“walkshed”. </p><br /> 356<p>In <a href="/posts/calculating-service-area-part-2-defining-the-walkshed">Part 2</a>, we<br /> 357defined what our walkshed is and defined how we would be linking it to the<br /> 358Census 2010. In this post, we’ll be doing the linking</p><br /> 359<p>I want to point out <a href="http://sproke.blogspot.com/2012/01/importing-2010-sf1-census-in-postgresql.html">sproke’s blogpost from<br /> 3602012</a><br /> 361where I got many of these ideas many years back.</p><br /> 362<h2>Parts</h2><br /> 363<ul><br /> 364<li><a href="/posts/calculating-service-area-part-1-introduction">Part 1: Introduction and Setup</a></li><br /> 365<li><a href="/posts/calculating-service-area-part-2-defining-the-walkshed">Part 2: Defining the Walkshed</a></li><br /> 366<li>Part 3: Demographics</li><br /> 367<li><a href="/posts/calculating-service-area-part-4-next-steps">Part 4: Next Steps</a></li><br /> 368</ul><br /> 369<h2>Importing the Census 2010 Data</h2><br /> 370<p>First we need to gather the data from <a href="https://www.census.gov/">census.gov</a> and<br /> 371initialized our database. I’m going to grab the documentation (sf1.pdf), along with<br /> 372some MS Access database files they use to store the schema.</p><br /> 373<pre><code>wget --continue -e robots=off --force-directories --no-parent https://www.census.gov/prod/cen2010/doc/sf1.pdf 374wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/04-Summary_File_1/SF1_Access2003.mdb 375wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/04-Summary_File_1/SF1_Access2007.accdb 376wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/04-Summary_File_1/Pennsylvania/ 377wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/05-Summary_File_2/SF2_MSAccess_2003.mdb 378wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/05-Summary_File_2/SF2_MSAccess_2007.accdb 379wget --recursive --continue -e robots=off --force-directories --no-parent https://www2.census.gov/census_2010/05-Summary_File_2/Pennsylvania/ 380</code></pre><br /> 381<p>Some notes, <code>--recursive</code> will grab linked files, <code>--continue</code> will try to<br /> 382start mid-file if it’s been partially downloaded. <code>-e robots=off</code> is naughty<br /> 383and tels wget to download things even if robots.txt tells you not to. (I don’t<br /> 384think it’s an issue here, I just add it all the time.) <code>--force-directories</code><br /> 385will create a directory structure begining with the domain so<br /> 386“https://example.com/a/b/c.d” will be downloaded to a file<br /> 387“example.com/a/b/c.d” — I find this makes keeping track of what things are and<br /> 388where they came from a lot easier. <code>--no-parent</code> tells wget not to go up the<br /> 389directory tree to get assets, which could cause it to begin going down sibling<br /> 390directories of the root we gave it.</p><br /> 391<p>Now, we are going import that data. The first step <br /> 392is to unzip the datafile. Be warned, the unziped data can swell in size. For<br /> 393example, the pa dataset is compressed 90%: the zip is 436MiB and uncompressed<br /> 394is 4.4GiB.</p><br /> 395<pre><code>unzip pa2010.sf1.zip 396</code></pre><br /> 397<p>Now that we have the files, let’s build the schema. Mostly because I love<br /> 398using the commandline whenever possible, I decided to stream the schema<br /> 399from the Access sample databases through some fix-ups and then into a file.</p><br /> 400<pre><code>mdb-schema SF1_Access2003.mdb | 401 # MS Access uses [] as idenifier quotes, let’s fix that 402 tr ‘[]’ ‘”“’ |</p> 403 404 # Some tables have a name suffix of `mod` because the first few have been 405 # removed to keep the table under 255 fields; let’s add them back in. 406 sed ‘/mod”$/{$!{N;s/mod”\n (/”\n (\n “<span class="caps">FILEID</span>” text,\n “<span class="caps">STUSAB</span>” text,\n “<span class="caps">CHARITER</span>” text,\n “<span class="caps">CIFSN</span>” text,/g}}’ | 407 408 # Remove spaces from identifiers. 409 sed -E ‘s/(\s+)”([<sup>”]+) ([</sup>”]+)”/\1\2_\3/g’ | 410 411 # Convert to PG type names 412 sed -E ‘s/Text\s*\([0-9]+\)/text/g’ | 413 sed -E ‘s/(Long )?Integer/bigint/g’ | 414 sed -E ‘s/Double/double precision/g’ | 415 416 # Table 45 is split into 2 parts; let’s merge them by deleting the trailing 417 # ); of the one and the <span class="caps">CREATE</span> <span class="caps">TABLE</span> along with the first few fields of the 418 # second. 419 sed ‘
419s/_PT1//g’ | 420 sed ‘/);$/{$!{N;N;N;N;N;N;N;N;s/.*_PT2.*LOG[^,]*//g}}’ | 421 422 # Just remove the quotes — they’re not needed anymore 423 tr -d ‘”’ > SF1_pg.schma.sql 424</code></pre> 425<p>Then, we can send that into postgres, but I want to wrap it all in a commit<br /> 426so that we don’t end up with a half-built database if something fails.</p><br /> 427<pre><code>cat &lt;(echo "set role census;") \ 428 &lt;(echo "begin; create schema census2010; set search_path = census2010;") \ 429 SF1_pg.schma.sql \ 430 &lt;(echo "commit;") | 431 psql -v ON_ERROR_STOP=1 census 432</code></pre><br /> 433<p>Now the fun part! We get to import the datafiles! I’m going to use the<br /> 434<a href="https://www.postgresql.org/docs/10/sql-copy.html">copy</a> statement and the <br /> 435psql <a href="https://www.postgresql.org/docs/10/app-psql.html#APP-PSQL-META-COMMANDS-COPY">\copy</a> meta-command. However, like above,<br /> 436I’m going to build up a sql file before running it so that everything can<br /> 437be wrapped in a transaction. Yes, this step can take a bit and if it fails<br /> 438towards the end it sucks, but in the end it doesn’t take that long (a few<br /> 439minutes) and it’s significantly easier than trying to scrub duplicates. (We<br /> 440haven’t built indices yet as having them during the input slows it down<br /> 441considerably.</p><br /> 442<pre><code>echo "begin;" &gt; import_script 443for input_file in pa000*2010.sf1; do 444 table_name=SF1_`echo $input_file | sed 's/pa\(.\{5\}\)2010.sf1/\1/'` 445 # Provides some feedback as `copy` will only output the number of rows imported. 446 echo "\\\echo $input_file $table_name" &gt;&gt; import_script 447 echo "\\\copy census2010.${table_name} from '${input_file}' with (format csv, header false);\n" &gt;&gt; import_script 448done 449echo "commit;" &gt;&gt; import_script
450cat import_script | psql -v ON_ERROR_STOP=1 census 451</code></pre><br /> 452<p>Now comes the more interesting part: the geoheader (links demographic records<br /> 453with the geography) is a fixed-width file (<em>i.e.</em> fields are not delimited, but<br /> 454are defined as between character positions on a line). What Sophia did was<br /> 455to build a staging table that would contain the raw lines, and then do a<br /> 456“insert into select” with a select statement that teased out the fields using<br /> 457the substring function. We’ll do a very similar thing. The fields are<br /> 458defined in the <a href="https://www.census.gov/prod/cen2010/doc/sf1.pdf">SF1</a><br /> 459documentation</p><br /> 460<pre><code>echo "set role census; CREATE TABLE census2010.geo_header_staging(data text);" | psql census 461echo "\\\copy census2010.geo_header_staging from 'pageo2010.sf1';" | psql census 462</code></pre><br /> 463<p>Then we’ll tease apart that staging table.</p><br /> 464<pre><code>INSERT INTO census2010.geo_header_sf1 (fileid, stusab, sumlev, geocomp, chariter, cifsn, logrecno, region, division, state, county, countycc, countysc, cousub, cousubcc, cousubsc, place, placecc, placesc, tract, blkgrp, block, iuc, concit, concitcc, concitsc, aianhh, aianhhfp, aianhhcc, aihhtli, aitsce, aits, aitscc, ttract, tblkgrp, anrc, anrccc, cbsa, cbsasc, metdiv, csa, necta, nectasc, nectadiv, cnecta, cbsapci, nectapci, ua, uasc, uatype, ur, cd, sldu, sldl, vtd, vtdi, reserve2, zcta5, submcd, submcdcc, sdelem, sdsec, sduni, arealand, areawatr, name, funcstat, gcuni, pop100, hu100, intptlat, intptlon, lsadc, partflag, reserve3, uga, statens, countyns, cousubns, placens, concitns, aianhhns, aitsns, anrcns, submcdns, cd113, cd114, cd115, sldu2, sldu3, sldu4, sldl2, sldl3, sldl4, aianhhsc, csasc, cnectasc, memi, nmemi, puma, reserved) 465SELECT 466nullif(trim(substring(data,1,6)), '') AS fileid, 467nullif(trim(substring(data,7,2)), '') AS stusab, 468nullif(trim(substring(data,9,3)), '') AS sumlev, 469nullif(trim(substring(data,12,2)), '') AS geocomp, 470nullif(trim(substring(data,14,3)), '') AS chariter, 471nullif(trim(substring(data,17,2)), '') AS cifsn, 472nullif(trim(substring(data,19,7)), '')::bigint AS logrecno, 473nullif(trim(substring(data,26,1)), '') AS region, 474nullif(trim(substring(data,27,1)), '') AS division, 475nullif(trim(substring(data,28,2)), '') AS state, 476nullif(trim(substring(data,30,3)), '') AS county, 477nullif(trim(substring(data,33,2)), '') AS countycc, 478nullif(trim(substring(data,35,2)), '') AS countysc, 479nullif(trim(substring(data,37,5)), '') AS cousub, 480nullif(trim(substring(data,42,2)), '') AS cousubcc, 481nullif(trim(substring(data,44,2)), '') AS cousubsc, 482nullif(trim(substring(data,46,5)), '') AS place, 483nullif(trim(substring(data,51,2)), '') AS placecc, 484nullif(trim(substring(data,53,2)), '') AS placesc, 485nullif(trim(substring(data,55,6)), '') AS tract, 486nullif(trim(substring(data,61,1)), '') AS blkgrp, 487nullif(trim(substring(data,62,4)), '') AS block, 488nullif(trim(substring(data,66,2)), '') AS iuc, 489nullif(trim(substring(data,68,5)), '') AS concit, 490nullif(trim(substring(data,73,2)), '') AS concitcc, 491nullif(trim(substring(data,75,2)), '') AS concitsc, 492nullif(trim(substring(data,77,4)), '') AS aianhh, 493nullif(trim(substring(data,81,5)), '') AS aianhhfp, 494nullif(trim(substring(data,86,2)), '') AS aianhhcc, 495nullif(trim(substring(data,88,1)), '') AS aihhtli, 496nullif(trim(substring(data,89,3)), '') AS aitsce, 497nullif(trim(substring(data,92,5)), '') AS aits, 498nullif(trim(substring(data,97,2)), '') AS aitscc, 499nullif(trim(substring(data,99,6)), '') AS ttract, 500nullif(trim(substring(data,105,1)), '') AS tblkgrp, 501nullif(trim(substring(data,106,5)), '') AS anrc, 502nullif(trim(substring(data,111,2)), '') AS anrccc, 503nullif(trim(substring(data,113,5)), '') AS cbsa, 504nullif(trim(substring(data,118,2)), '') AS cbsasc, 505nullif(trim(substring(data,120,5)), '') AS metdiv, 506nullif(trim(substring(data,125,3)), '') AS csa, 507nullif(trim(substring(data,128,5)), '') AS necta, 508nullif(trim(substring(data,133,2)), '') AS nectasc, 509nullif(trim(substring(data,135,5)), '') AS nectadiv, 510nullif(trim(substring(data,140,3)), '') AS cnecta, 511nullif(trim(substring(data,143,1)), '') AS cbsapci, 512nullif(trim(substring(data,144,1)), '') AS nectapci, 513nullif(trim(substring(data,145,5)), '') AS ua, 514nullif(trim(substring(data,150,2)), '') AS uasc, 515nullif(trim(substring(data,152,1)), '') AS uatype, 516nullif(trim(substring(data,153,1)), '') AS ur, 517nullif(trim(substring(data,154,2)), '') AS cd, 518nullif(trim(substring(data,156,3)), '') AS sldu, 519nullif(trim(substring(data,159,3)), '') AS sldl, 520nullif(trim(substring(data,162,6)), '') AS vtd, 521nullif(trim(substring(data,168,1)), '') AS vtdi, 522nullif(trim(substring(data,169,3)), '') AS reserve2, 523nullif(trim(substring(data,172,5)), '') AS zcta5, 524nullif(trim(substring(data,177,5)), '') AS submcd, 525nullif(trim(substring(data,182,2)), '') AS submcdcc, 526nullif(trim(substring(data,184,5)), '') AS sdelem, 527nullif(trim(substring(data,189,5)), '') AS sdsec, 528nullif(trim(substring(data,194,5)), '') AS sduni, 529nullif(trim(substring(data,199,14)), '')::double precision AS arealand, 530nullif(trim(substring(data,213,14)), '')::double precision AS areawatr, 531nullif(trim(substring(data,227,90)), '') AS name, 532nullif(trim(substring(data,317,1)), '') AS funcstat, 533nullif(trim(substring(data,318,1)), '') AS gcuni, 534nullif(trim(substring(data,319,9)), '')::bigint AS pop100, 535nullif(trim(substring(data,328,9)), '')::bigint AS hu100, 536nullif(trim(substring(data,337,11)), '') AS intptlat, 537nullif(trim(substring(data,348,12)), '') AS intptlon, 538nullif(trim(substring(data,360,2)), '') AS lsadc, 539nullif(trim(substring(data,362,1)), '') AS partflag, 540nullif(trim(substring(data,363,6)), '') AS reserve3, 541nullif(trim(substring(data,369,5)), '') AS uga, 542nullif(trim(substring(data,374,8)), '') AS statens, 543nullif(trim(substring(data,382,8)), '') AS countyns, 544nullif(trim(substring(data,390,8)), '') AS cousubns, 545nullif(trim(substring(data,398,8)), '') AS placens, 546nullif(trim(substring(data,406,8)), '') AS concitns, 547nullif(trim(substring(data,414,8)), '') AS aianhhns, 548nullif(trim(substring(data,422,8)), '') AS aitsns, 549nullif(trim(substring(data,430,8)), '') AS anrcns, 550nullif(trim(substring(data,438,8)), '') AS submcdns, 551nullif(trim(substring(data,446,2)), '') AS cd113, 552nullif(trim(substring(data,448,2)), '') AS cd114, 553nullif(trim(substring(data,450,2)), '') AS cd115, 554nullif(trim(substring(data,452,3)), '') AS sldu2, 555nullif(trim(substring(data,455,3)), '') AS sldu3, 556nullif(trim(substring(data,458,3)), '') AS sldu4, 557nullif(trim(substring(data,461,3)), '') AS sldl2, 558nullif(trim(substring(data,464,3)), '') AS sldl3, 559nullif(trim(substring(data,467,3)), '') AS sldl4, 560nullif(trim(substring(data,470,2)), '') AS aianhhsc, 561nullif(trim(substring(data,472,2)), '') AS csasc, 562nullif(trim(substring(data,474,2)), '') AS cnectasc, 563nullif(trim(substring(data,476,1)), '') AS memi, 564nullif(trim(substring(data,477,1)), '') AS nmemi, 565nullif(trim(substring(data,478,5)), '') AS puma, 566nullif(trim(substring(data,483,18)), '') AS reserved 567FROM census2010.geo_header_staging; 568</code></pre><br /> 569<p>Now, let’s delete that staging table.</p><br /> 570<pre><code>drop table census2010.geo_header_staging; 571</code></pre><br /> 572<p>Let’s build some indices and foreign keys.</p> 573<pre><code>echo “begin;” > create_index 574for table in `echo “\\\dt census2010.*” | psql -tA census | grep -v data | cut -f 2 -d \|`; do 575 cat <<<span class="caps">EOS</span> >> create_index 576create index on census2010.${table}(fileid); 577create index on census2010.${table}(stusab); 578create index on census2010.${table}(chariter); 579create index on census2010.${table}(cifsn); 580create unique index on census2010.${table}(logrecno); 581 582 <p>alter table census2010.${table} alter column fileid set not null;<br /> 583alter table census2010.${table} alter column stusab set not null;<br /> 584alter table census2010.${table} alter column chariter set not null;<br /> 585alter table census2010.${table} alter column cifsn set not null;<br /> 586alter table census2010.${table} alter column logrecno set not null;</p> 587 588 <p>alter table census2010.${table} add primary key using index ${table}_logrecno_idx;</p> 589 590 <p>EOS<br /> 591done<br /> 592echo “commit;” >> create_index<br />
593cat create_index | psql -v ON_ERROR_STOP=1 census</p> 594 595 <p>create unique index on census2010.geo_header_sf1(fileid, stusab, chariter, cifsn, logrecno);</p> 596 597 <p>echo “begin;” > create_fk<br /> 598for table in `echo “\\\dt census2010.*” | psql -tA census | grep -v data | grep -v geo | cut -f 2 -d \|`; do 599 cat <<<span class="caps">EOS</span> >> create_fk<br /> 600alter table census2010.${table} add constraint ${table}_geo_fk foreign key (logrecno) references census2010.geo_header_sf1(logrecno);<br /> 601EOS<br /> 602done<br /> 603echo “commit;” >> create_fk<br /> 604cat create_fk | psql -v ON_ERROR_STOP=1 census<br /> 605</code></pre><br /> 606<h2>Service Area Population</h2><br /> 607<p>In order to speed up lookups between the <span class="caps">TIGER</span> data, let’s create an index to<br /> 608use. <code>sumlev</code> (Summary Level) is the aggregation that the row represents,<br /> 609<em>e.g.</em> state, county, track, block group, block, &c. They’re defined in the<br /> 610<a href="https://www.census.gov/prod/cen2010/doc/sf1.pdf">SF1</a> documentation. <code>101</code> is<br /> 611the block level.</p><br /> 612<pre><code>create index on census2010.geo_header_sf1 (state, county, tract, block) where sumlev = '101'; 613</code></pre><br /> 614<p>Calculating the total population of the county can be done by simply using<br /> 615field <code>p0010001</code> (total population) in table <code>sf1_00001</code> (also defined in the<br /> 616SF1 documentation), but only for the rows representing the block summary level.<br /> 617Without filtering on the summary level the population would be counted for<br /> 618every summary level in the county.</p><br /> 619<pre><code>select sum(p0010001) 620from census2010.sf1_00001 621join census2010.geo_header_sf1 using (logrecno) 622where county = '003' 623 and sumlev = '101'; 624</code></pre><br /> 625<p>For the population in our walkshed, we will use the <code>stops_walkshed_blocks</code><br /> 626table that we generated in Part 2. As a recap, that is every block that is<br /> 627adjacent to a segment of road that is within a quarter mile of a bus stop. (The<br /> 628reason for the <code>group by</code> is that <code>stops_walkshed_blocks</code> is broken down by<br /> 629stop, which obviously leads to some blocks being counted multiple times.)</p><br /> 630<pre><code>select sum(p0010001) 631from ( 632 select max(p0010001) as p0010001 633 from census2010.sf1_00001 634 join census2010.geo_header_sf1 using (logrecno) 635 join tiger2017.tabblock 636 on statefp10 = state 637 and countyfp10 = county 638 and tractce10 = tract 639 and blockce10 = block 640 join stops_walkshed_blocks 641 on stops_walkshed_blocks.geoid10 = tabblock.geoid10 642 where county = '003' and sumlev = '101' 643 group by tabblock.geoid10) x; 644</code></pre><br /> 645<h2>Results</h2><br /> 646<p>This yields a total population of 1,223,348 and a service area population of<br /> 647775,835 which is just over 63% of the county’s population.</p></p> 648 649 </div> 650 </div> 651 652 <hr> 653 <footer class="footer" id="footer"> 654 <div 655 <p>This site is made possible by <a href="https://github.com/jimktrains/gus">gus</a> and <a href="http://twitter.github.com/bootstrap/">Bootstrap</a>. <a href="http://jimkeener.com/rss.xml">RSS</a></p> 656<p>My words are copyright Jim Keener 2023. The technology is all <a href="http://www.gnu.org/philosophy/free-software-for-freedom.html">F/OSS</a>.</p> 657 </footer> 658 </div> <!-- /container --> 659 </body> 660</html>
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.