Showing posts with label QGIS. Show all posts
Showing posts with label QGIS. Show all posts

Sunday, August 9, 2026

QGIS 4.0 to 4.2.0 to 4.2.1 in July heat with NetBSD 11.0 hot off the digital press

 Around the first of Jul 2026 I was building QGIS 3 or 4 on 3 NetBSD systems, 1 fast, 2 slow, in order to check that each worked with PostgreSQL 14, 16 or 18, respectively. It mostly worked, albeit as slow as molasses in winters of yore.

To do the builds requires pkgsrc source code, as the pre-packaged versions are only built with version 16 database driver/library support.

CPU temperature + Processor load 24 hours+

On the 6th, I reworked the Bengies drive-in geocache and tent map.


The Bengies Drive-In

On the bigger system, 4.0 QGIS

QGIS 4 splash screen

Next thing I know, 4.2.0 is sitting in work-in-progress to I backed up and restarted the build.





Very nice, and worked about the same as the 4.0.

4.2.0 About Screen


I tried the Authentication option briefly, didn't get screen shots, and then. poof again, 2026Q2 pkgsrc became available after bulk builds, NetBSD 11.0 was released after a BETA and 6 or more Release Candidates.


By 31-July, a geospatial PDF:






When I finished 11.0 upgrades and wiped pkgsrc for a clean slate, there came one more work-in-progress version, 4.2.1.

Or "final".

QGIS 4.2.1 on NetBSD 11.0



Monday, July 20, 2026

QGIS 3x3 on NetBSD 10/11 + PostgreSQL 14/16/18

 

Started making QGIS version 3.44.8 early June 2026 on NetBSD processors

jump ahead:

$ ps auxww | head -2

USER        PID %CPU %MEM      VSZ    RSS TTY    STAT STARTED     TIME COMMAND

root      24310 93.3  1.1   427668 355252 pts/6  O+   12:18AM  0:02.14 /usr/pkg/gcc12/libexec/gcc/x86_64--netbsd/12.4.0/cc1plus -quiet -I /usr/pkgsrc/wip/qgis/work/QGIS-final-4_0_3/cmake-pkgsrc-build/src/core/qgis_core_autogen/include

Wed Jun 10 00:25:37 UTC 2026

$ TZ=EST date

Tue Jun  9 19:25:45 EST 2026

=> Checking file-check results for qgis-3.44.8
=> Creating binary package /usr/pkgsrc/geography/qgis/work/.packages/qgis-3.44.8.tgz


/// pkg_delete postgresql18-client qgis-3.44.8  py313-psycopg2-2.9.11


$ pkgin  show-deps "qgis" | grep -i client
        postgresql16-client>=16.0nb1


https://pkgsrc.se/wip/qgis

Version history: (Expand)

    (2026-06-02) Updated to version: qgis-4.0.3
    (2026-05-10) Updated to version: qgis-4.0.2.18
    (2026-03-09) Updated to version: qgis-4.0.0.0
    (2026-02-22) Updated to version: qgis-3.99.3
    (2026-02-21) Updated to version: qgis-3.99.2

FREEBSD ~~~

foobar:/etc # pkg search qgis
qgis-3.44.11                   Free and Open Source Geographic Information System
qgis-ltr-3.40.15_1             Free and Open Source Geographic Information System

something bollixed. Back to NetBSD.
~~~

amd64-12-core
$ ls -ltr /usr/pkgsrc/packages/All/qg*
-rw-r--r--  1 root  wheel  123201876 Jun  9 17:15 /usr/pkgsrc/packages/All/qgis-3.44.8.tgz
-rw-r--r--  2 root  wheel  124614138 Jun 26 16:10 /usr/pkgsrc/packages/All/qgis-3.44.11.tgz
-rw-r--r--  2 root  wheel  127333035 Jul 13 11:34 /usr/pkgsrc/packages/All/qgis-4.0.3.tgz
-rw-r--r--  2 root  wheel  131910253 Jul 15 23:27 /usr/pkgsrc/packages/All/qgis-4.2.0.tgz

amd64-4-core
5.3$ ls -l /root/pkg-local/qg*
-rw-r--r--  1 root  wheel  127333035 Jul 13 20:02 /root/pkg-local/qgis-4.0.3.tgz

amd64-2-core
[     1.000004] cpu0 at mainbus0 apid 0
[     1.000004] cpu0: Use lfence to serialize rdtsc
[     1.000004] cpu0: Intel(R) Core(TM)2 Duo CPU     L9400  @ 1.86GHz, id 0x1067a
[     1.000004] cpu0: node 0, package 0, core 0, smt 0
[     1.000004] cpu1 at mainbus0 apid 1
[     1.000004] cpu1: Intel(R) Core(TM)2 Duo CPU     L9400  @ 1.86GHz, id 0x1067a
[     1.000004] cpu1: node 0, package 0, core 1, smt 0
[     1.000004] acpi0 at mainbus0: Intel ACPICA 20221020



2 core had faults, see:

RuntimeError: NumPy was built with baseline optimizations:
(X86_V2) but your machine doesn't support:
(X86_V2).

Pictures and exhibitions

Kept track of long compile times (~18 hours) on the slowest machine, as well as faster, looking at CPU temperature and process load (distributed across 1 or many "cores").

[
mastodon cc me
https://assets.chaos.social/media_attachments/files/116/846/959/384/561/931/original/96c70f8d89d6e49f.png

Go Speed Racer Go!


Pre 3-way build side quest:


Initial launch 4.0.3 June 2026


-rw-r--r--  1 jim  jim 20978 Jun 10 12:46 GIS/spots-rooms-q4.qgz



Side quest:

QGIS 3.44 on the machine I chose did not build, so the check with PG SQL 14 was inadequate. I changed the SQL version on that machine to 18, and it runs with the broken CPU message disabling various plugins (yet to be found).

This is the override of the database version when building from source, not installing pre-fabs.

$ pwd
/usr/pkgsrc/mk

$ diff pg*
28,30c28
< # PGSQL_VERSION_DEFAULT?=             16
< # Fri Jun 26 04:32:43 UTC 2026
< PGSQL_VERSION_DEFAULT?=               18
---
> PGSQL_VERSION_DEFAULT?=               16


Screengrab needs to be built from scratch; what I previously did with an X capture app has fallen into disrepair. At least xli still works.

$ ls -l `which qgis`
-rwxr-xr-x  1 root  wheel     390320 Jul  3 12:43 /usr/pkg/bin/qgis

-rw-r--r--  1 root  wheel  124615383 Jul  3 12:45 ./pkg-local/qgis-3.44.11.tgz

$ psql --version
psql (PostgreSQL) 18.4


Post 3-way build side quest of 3x3, going to 4:

Back to the fast machine and the newest (l+g). The work-in-progress version went from 4.0.3 to 4.2.0. Wow, quick change!

Build time ~8 hours with or without fumbles. Needs to be redone with a clean slate and only 2026Q4 pkgsrc to get things level. But it starts and connects to databases as it should.

$ time qgis --version
QGIS 4.2.0-Belém do Pará 'Belém do Pará' (exported)

real    0m2.522s
user    0m2.458s
sys     0m0.039s







Quick duck back to 3.44.11



Comparison table of the 12-core with 4.2.0 and the 2-core with 3.44.11.



Sidebar issue with Python caused the active plugins to be empty on the 2-core, but otherwise they are pretty much the same beyond QGIS itself and the Qt version going from 5 to 6.

Neither system shows a value for SFCGAL or GeographicLib. Not sure what impact that has.

Working on a list of projects for QGIS 4.2, including waypoints of Historic Trails, paths in the woods, and geocaching at the Bengies Drive-In Theatre.

Monday, June 30, 2025

Take the 40 Bus Route Out Route 40

 After getting feedback, advice, and counsel about my recent public transport mapping updates in OpenStreetMaps, I started adding the missing "40" bus route to the Baltimore Metro Area. Part of the planning was whether to include bus stops and stopping locations, part was identifying errors and omissions in the OSM base, and part was being able to visualize current routes with OSM data.


Then a fragment:



I can see the 40 on existing bus stop signs:




This stop includes the OR (Orange), the 40, 59, 62, and 160 routes. When I looked, the OSM data for stop 7597 showed '59;62;160;OR' while the MTA showed '59;62;160;OR', as if the bus went to the Essex Park + Ride on the other side of the street instead of this stop.


So, I am unsure what the right data should be, until I ride the bus, or observe at the stop.

The guidance for adding a new route ("Step 1 - Make sure that each bus stop in the route has been added to OpenStreetMap") said check all the nodes. Given the above discrepancy, I hope to add at least one-way next, as small fixes could happen after, like when the routes/schedules change.

The last overview map update page from MTA [ mta.maryland.gov/transit-maps ] is from 2022, and the 40 route is not on the "BaltimoreLink Interactive System Map". That's a map from Google rather than hosted by the government. However, extrapolating shows this page is available in 2025:


Methodology

For myself, and others who may wish to get some GIS mojo working, general steps I've taken to gather data, analyze it, and visualize in various ways to find improvements to OSM content. The earlier post [ https://jspath55.blogspot.com/2025/06/which-bus-such-bus.html ] covers some of this; apologies for repeating myself.

MTA Data

The Maryland Department of Transportation includes the Mass Transit Administration. In addition to the schedule site listed above, there is a GIS portal with bus stops and routes, though I've not gotten much out of the latter.


Home > services > Transportation > MD_Transit (FeatureServer)


Then:


And end up here:



GIS ARC

Part 2 (right side)


The way to pull the data is via:

https://geodata.md.gov/imap/rest/services/Transportation/MD_Transit/FeatureServer/9/query

https://geodata.md.gov/imap/rest/services/Transportation/MD_Transit/FeatureServer/9/query?where=stop_id+is+NOT+NULL ...


"All" query: stop_id is NOT NULL

Output fields: stop_id, stop_name, Routes_Served, shelter

Results in a linear format:

# records: 4549

stop_name: Cylburn Ave & Greenspring Ave
Routes_Served: 31, 91, 94
Shelter: Yes
stop_id: 1
Point:
X: -8533795.9128
Y: 4772066.8330999985

stop_name: Lanier Ave & Sinai Hospital
Routes_Served: 31, 91, 94
Shelter: No
stop_id: 3
Point:
X: -8534126.0864
Y: 4772153.208300002

[...]

I parsed this into a CSV format to load into QGIS. The wrinkle was finding how to deal with the X Y coordinates returned rather than the expected latitude/longitude pairs. QGIS was happy when I set the coordinate reference system (CRS) to "ESRI:102100 - WGS_1984_Web_Mercator_Auxiliary_Sphere."

OSM data

I used OverPass-Turbo.eu to get extracts of bus stopping locations and bus stops, then separate them.

[out:csv("ref", "public_transport", "name", "route_ref", "highway", "network", "operator", ::lat, ::lon, ::id)];

(
  {{geocodeArea:"Maryland"}};

  )->.metroArea;

node(area)[bus=yes];
out;

I used "Maryland" rather than "Baltimore" to pick up Anne Arundel County and other nearby locales, though this expanded the return to include eastern and western counties as well as DC Metro. Better too much data than not enough here.

The first tries did not include position or OSM node ID; eventually I found the special names needed to get theses critical columns. After splitting into 2 files to allow QGIS layers for each, over 4,000 bus stops found though under 1,000 stopping locations.

$ wc -l OSM_platforms_20250623.csv OSM_stops_20250623.csv
  4240 OSM_platforms_20250623.csv
   928 OSM_stops_20250623.csv
  5168 total

Sample export:



Add these 2 layers, plus the MTA extract, and start the reconciliation.

Things found:

  • MTA and OSM sequence the routes differently, and the MTA records appear in yet another order compared to bus stop signs. I wrote a Perl script to convert the MTA data into the order OSM expects to make it easier to update wrong data.
  • The bus stop names differ in upper case/lower case, and various abbreviations ("nb" versus "northbound", for example.
  • Typos in both OSM and MTA data conglomerating 2 or more routes to one number ("78150" instead of "78, 150" or "78;150".
  • Null data nodes: apparent OSM bus stops without details (e.g. openstreetmap.org/node/2743861086 ). These are presumably orphans, though they may be incomplete nodes that should remain. If no MTA stops nearby, fix me candidates.
  • Typos of other sorts, missing stops, incomplete or duplicate descriptions.

Things fixed:

I alternated trying to simplify the updates to be made by only looking at platforms, but decided that version 2 with a stopping location paired with the platform is the preferred way. I created a spreadsheet to have a local index of the 40 bus routes in either direction. The columns included nodes for each type, stop numbers and names, and an alternate name per the official signage ("Eastpoint" as the area for Eastern & 54th, east or west bound).

Only a handful of bus stop areas include a stopping location, so that is my current task. After that, on to steps 2 and 3:

Step 2 - Create the new bus relations
Step 3 - Create a route master relation


The QGIS map view, overlaying MTA, OSM stops/platforms:



By migrating these layers into a PostGIS/PostgreSQL database (using a CRS under 10,000), I can now write SQL queries with outer joins and other comparison techniques.

Stops in MTA but not OSM:

select
  ref,
  substr(trim(route_ref),1,8) as routes,
  trim(name) as name
from
  "OSM_platforms_20250623"
where
  cast(ref as text) NOT in (
    select
      field_4
    from
      "bus-stops-20250611"
  )
;
[...]
   14166 | 30;34    | Northern Parkway & Park Heights Avenue Westbound Far-side
   14215 |          | Tradepoint Atlantic Floor & Decor
   14216 |          | Tradepoint Atlantic Royal Farms
  190036 |          | Baltimore Greyhound Terminal
 3002008 |          | Good Luck Rd at Oakland Ave
 3003202 |          | Good Luck Rd at Parkdale High School
 3003205 |          | Good Luck Rd at Oakland Ave
(133 rows)


Stops in MTA but not OSM:

select
      upper(substr(mta.field_1,1,40)) as nm,
      substr(mta.field_2,1,32) as f2,
      mta.field_4,
      osm.ref
from
      "bus-stops-20250611" mta
left outer join
      "osm_platforms_20250623_md" osm on mta.field_4 = cast(osm.ref as text)
where
  osm.ref is NULL
or
  cast(osm.ref as float) <= 0
order by
  cast(mta.field_4 as float)
;

This one included DC Metro areas, so needs more tuning to be useful.

Show mismatched route lists per stop:

select
      upper(substr(mta.name,1,40)) as nm,
      substr(mta.route_osm_style,1,32) as mta_route_osm,
      substr(osm.route_ref,1,32) as osm_route,
      mta.stop,
      osm.ref
from
      "bus-stops-20250611+0626" mta
left outer join
      "osm_platforms_20250623_md" osm on cast(mta.stop as text)  = osm.ref
where
  (osm.ref is NULL or cast(osm.ref as text) >= '0')
and
  mta.route_osm_style != osm.route_ref
order by
  mta.stop
;

The later version of the MTA list includes the reordered route list, as comparing the strings in SQL was not an obvious approach.  Example results:

[...]
      upper(substr(mta.name,1,40)) as nm,
      substr(mta.route_osm_style,1,32) as mta_route_osm,
      substr(osm.route_ref,1,32) as osm_route,
      mta.stop
[...]
LOCH RAVEN BLVD & STATE HWY DR SB 53 53;103;104 422
LOCH RAVEN BLVD & SAYWARD AVE SB 53;GR 53;103;104;GR 435
SAINT PAUL ST & 31ST ST SB 51;95;SV 95;SV 482
SAINT PAUL ST & 29TH ST SB 51;95;SV 95;SV 483
CHARLES ST & REDWOOD ST NB 51;65;67;95;103;410;411;420;GR;S 51;65;67;95;103;164;410;411;420; 506
SAINT PAUL ST & FAYETTE ST FS SB 67;76;95;103;120;210;215;310;410 67;76;95;103;120;164;210;215;310 516
CHARLES ST & HYATT NB 51;67;94;95;103;SV 51;67;94;95;103;103;164;SV 518
CHARLES ST & FAYETTE ST NB 51;95;103;410;411;420;GR;SV 51;95;103;410;411;420;425;GR;SV 519
CHARLES ST & READ ST NB 51;95;103;GR;SV 51;95;103;410;GR;SV 525

I photographed the bus stop sign at Charles and Read Streets (northbound) only yesterday, confirming the list MTA had. And updated the node: openstreetmap.org/node/6254842377

Charles + Read Streets Baltimore

U.S. Route 40, as opposed to MTA Bus Route 40, is the "National Pike"; see also wikipedia.org/wiki/Old_National_Pike.



Tuesday, June 10, 2025

Which Bus? Such bus.

 

It started innocently when I noticed a local transit shelter had been moved to the other side of the highway but that wasn't showing up in Open Street Maps. When I looked deeper, I found OSM had 1 shelter and 3 bus stop like nodes where the bus formerly stopped.

Because I was inexperienced at using OSM for public transport routing as opposed to vehicle directions, I knew I should be careful making edits less valid data or relationships be lost or confusing to others.

It later dawned on me that this wasn't 1 stop being removed, or moved, it was 2 stops replaced by one nearby. So I chose to delete the 3 bogons, and move the leftover down the road. I think I did it right and am waiting for various map sites to get these updates and see how they look.

But like pulling a loose thread that turns into an unravel-fest, I found more "errors and omissions" on other stops beyond where I started. Primary bugs:

  • Same bus stop ID number on 2 or more nodes. Should not happen.
  • Some stops had a route identifier; some did not.
  • Some stops had links/routes showing line numbers/names; some did not.
  • Nodes were classified as "bus" "bus stop" "bus station" willy-nilly
  • Multiple routes at one stop were mis-keyed as unseparated values
  • A new route is not on the map; does it exist?


Whispering Woods transit stop, in the area now known as Hopkins Point.



Bus stops seem duped.

Insert oops here

@neatnit@fosstodon.org Thank you for the guidance. Seems redundant stops were added by different people over time, and I was uncertain best approach. This stop was moved to the other side of the highway, so passengers go east one stop to a loop, then go west.


[out:json][timeout:25];
// gather results

/*
*/

node["highway"="bus_stop"]["ref"]({{bbox}});
for(t["ref"]) {
  if(count(nodes) > 1) {
    out tags center;
  }
}

// print results
out geom;


fix check:

https://overpass-api.de/api/sketch-line?network=BaltimoreLink&ref=59&operator=



The above is one view of the metro bus routes, showing the phantom (new) route #40.

It's ARCGIS, with Maryland State as the client/portal.

Another view:


Also ESRI-based, but with a more useful diagram of the route *and* stops. And termini.

https://www.mta.maryland.gov/schedule/stops/40

The MTA has a not-quite-current list of all bus stops, from a non-GUI menu path:

https://geodata.md.gov/imap/rest/services/Transportation/MD_Transit/FeatureServer/9/query

Simple-enough query-by-example fill-in-the boxes:


SELECT
stop_id, stop_name, Routes_Served, shelter
WHERE 
stop_id is NOT NULL

which returns columnar date in separate rows, plus coordinates.

e.g.

stop_name: OLD EASTERN AVE & SELIG AVE fs eb
Routes_Served: OR, 59, 62
Shelter: No
stop_id: 4522
Point:
X: -8511670.83
Y: 4766621.937700003

stop_name: OLD EASTERN AVE & ESSEX AVE eb
Routes_Served: OR, 59, 62, 160
Shelter: No
stop_id: 4524
Point:
X: -8511141.9511
Y: 4766737.762000002

(raw:)
/stop_name: / EASTERN BLVD & VIRGINIA AVE eb
/Routes_Served: / OR, 59, 62, 160
/Shelter: / No
/stop_id: / 7597
/Point:/
/X: / -8514337.1545
/Y: / 4765238.0484

The coordinate reference system I found worked with QGIS is "ESRI meters", also known as:


PROJCRS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
    ID["ESRI",102100]]

I used the findings of the duplicate stop search to check location and routes, as well as naming standards (which someone took care of 99.4% not too many years ago).






The zip-ties are an indicator of datum freshness. Broken signboard likewise.

Next steps:

  • Add the missing 40 bus to OpenStreetMaps as a blank spot on the canvas.
  • Find duplicate stops by checking how close together they are, that aren't headed in opposite directions.
  • Figure out what to do about the majority of missing routes (as relations) for nodes.
  • Other errors and omissions as they appear, i.e. quality assurance.
  • Share timetable/schedule/interval info. See: chaos.social/@jspath55/114659131418285869
  • Compare recorded comfort features (benches/shelters) with ground conditions.








Wednesday, April 9, 2025

Food Service Mapping

Food banks run by charitable organizations are a public service where private entities fill in where government runs short. Picking up, storing, and sharing food throws lifelines to those in need. I volunteered to create updated maps for Scouting America and decided to use QGis and related software tools.

(1) OSM

Going from the lowest level up, I added OpenStreetMaps (OSM). This is an easy drag-and-drop from the base QGIS sources into the current project. Depending on the desired result, I change the transparency from none to 50% more or less.

The standard OSM layer sources from openstreetmap.org. I've found, using tools like Viking, that variations on the main source can be used, and I prefer to try the Humanitarian one. Eventually I found the wiki page which let me set the feed, like this:

https://wiki.openstreetmap.org/wiki/Raster_tile_providers

https://a.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png

(2) Org

The second layer is the organization table. These have street addresses in the source, so I geo-located them with a Census Bureau tool that works via command line.

Once located, the QGIS data are converted to geometry. In a PostgreSQL database, I imported organization details, which had only street addresses not geo-data. To convert from the address to a lat/lon point I found a Census page that I could script.


$ lynx -dump "https://geocoding.geo.census.gov/geocoder/locations/address?street=10001%20Bird%20River%20Rd&state=md&zip=21220&benchmark=2020"  | grep Interpolate
   Interpolated Longitude (X) Coordinates: -76.432431362395
   Interpolated Latitude (Y) Coordinates: 39.356117523642

Ran the conversion steps and viewed points and polygons.  Once the latitude and longitude are set, the PostGIS function to make this into valid point(s) I used is:

UPDATE org.org SET geom = ST_SetSRID(ST_MakePoint(lng, lat), 4326);

(3) Tract shapes

Next level layer shows US Census tracts. I found a couple sources of these shape files, and used the Maryland State data.

https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2020&layergroup=Census+Tracts


[omitting how to add this layer from ZIP or shape files]

To match the boundaries of the local district I wrote a filter, first including only Baltimore County tracts and next excluding specific tracts not in scope.

QGIS lets me save and load a filter definition, more or less a SQL code snippet.

<Query>"COUNTYFP"='005' and "tractce"  NOT IN (
        '400100',
        '400200',
        '400400',
[...]
        '494201',
        '494202',
        '980000',
        '980100',
        '980200',

        '9999999999'
)
</Query>

The blank line prior to the query end lets me run a unique sort to more quickly find gaps and add or subtract tracts.

(4) Unit Tracts Coverage

The coverage mapping uses a simple table with tract and unit/organizations. In order to connect these to the tracts I created a view.

Thursday, March 21, 2024

WHERE clause, or FILTER? QGis and small databases

Continuing with getting more QGIS skills, I planned to run the tutorial for changing symbols based on rules or other data conditions. Before I got there, I found a treasure trove of GIS data feeds from the local state government, particularly apt as the sources include property shapes as well as many details on property values. I had looked at the Comptroller's site first only to find the better results on the state planning site.

I had gathered a few data records on specific properties from public information about drinking water wells and mapped a good set of views with QGIS. I was trying not to trace property lines from downloaded imagery and was rewarded for my delay by finding the official records. Of course with caveats, as with any data feed, especially government fed.

I found 4 data sets eventually, 3 as data points and the fourth as polygons. 2 of the 3 point files have  "property" records with one having few columns and the other many; the 3rd has building-specific info.

For my first data loads, I unzipped the procured files then added that data set as a shapefile by telling QGis that the file was the one with the dBASE extension (.dbf). Later I tried point the "Add Layer" file finder at the ZIP file without unpacking the archive and that seemed to work fine. Yay QGis!



Having the data in the GIS console opened many more possibilities. Once I figured out which key column to use (Account ID, go figure) I was able to create a join in QGis. A different experience than mine with database joins with command line SQL scripts, and a bit wrinkly around the edges I'd say.

Part of what I learned in the loading attempts was that source file names might need to be escaped in order to write SQL scripts or filters with QGis itself. The couple files with ASCII/mixed case simply worked with quotes around it. The file name with a Unicode "em dash" on the other hand i struggled with.

SELECT * FROM "public"."CAROparcels0124 — CAROMDPV.shp" LIMIT 10

I worked on creating a view off that table only to hit syntax errors due to bad encoding. I had backslash feedback of multiple kinds "\x0ef" e.g. Eventually I used a connection that maintained the correct characters to record the view definition.

The coolest feature that I found in QGis in this area was being able to drag a locally sourced shape file matrix layer with the GUI and drop it on top of an already open database connection. And hey! The DDL (data definition language) is complete. No scripting, no column defining, just done. I tried to do a screen grab while the create table and insert statements were running, but even with a Raspberry Pi as the "database server."


This showed the stage at 93% complete. I learned to be patient and just let the gears run.

I went back to database filters quickly after checking out QGis filter menu, for performance reasons. If I was looking at a wide area I would want to be able to pan from one small space to another at the price of slow start-up and possibly slow pans as more data are pulled into active stores. On a quick machine you may not notice such sluggishness; for learning purposes I am trying the tutorial lessons on high-ish end x86s and on low-ish end Raspberry Pis. 15,000 records transferred with 100 columns per row is noticeably a lag.

Knowing the data helps decide which fields are filter candidates. In this effort, I recognized some of the columns by their content, not so much their names, and I'll eventually read the data dictionary docs to become more familiar. The interesting ones included kinds of industrial and housing/zoning classification, year structure built and whether a building was on site.

I ended up choosing a political jurisdiction code after trying census tract and looking at other candidates. The first 2 characters are the county, which is constant in the files I pulled, with the next 2 being a district (state level?). That reduced the record counts from 15,000 or so to 3,000, depending on which source file I filtered.

The view definition gets all columns (for simplicity); if I had to, I would look for columns with all NULL or single-values and omit them.

create or replace view carobldg_dist07_v as
select * from "CaroBldg"
where acctid like '0607%'
;

The joins seemed pretty straightforward until I found one property with 2 IDs. Or one ID and 2 property uses. I forget now, and this is my *not to self* to run Pareto over there.

In a later run, the Baltimore County file has 10x what Caroline has, meaning 1/10 the speed in general. Again, character encoding issues.


(~4 minutes for one copy/paste)

> select count(*)  from "BACOparcels0124 \U+00E2\U+0080\U+0094 BACOPOLY.shp"
;
 count  
--------
 233227
(1 row)




2024-03-16T01:02:21     INFO    Task complete : Exporting BACOparcels0124 — BACOMDPV.shp2024-03-16T01:06:08     INFO    Task complete : Exporting BACOparcels0124 — BACOPOLY.shp


And, below, coordinate reference misses; more on that later once I grok the state-provided X-Y pairs.

Below is the Pi4 image of QGis 3.28 on NetBSD 10.0 (RC6).
The log times don't show exactly as I'd like; still looking for any performance logs.
I think that message about shapes is for the data load step, except they aren't in order it seems.

The highlighted area with data details on the right side is part of Baltimore's drinking water supply, showing a border of the City-owned property. Between that 500+ acre parcel is a thin green strip labelled "Loch Raven Reservoir Cooperative Wildlife Management Area."
On the left, I now have a stack of data sources, including "all" of the county info, or a portion via a database view. I'd plan on using the large source, trying different filter logic, and creating more views for the useful fractions.


(Santa, please bring us geoPDF in a future pkgsrc qgis package!)

Monday, March 4, 2024

Screen Scraping GIS Tax Ditch Lines Substandard

Better to get the GIS data someone else has already digitized than to try to digitize from a screenshot. That is what I tried a few months back before I found better source data recently delineating where stormwater drainage goes.

The idea that ditches are taxable, or tax exempt, is still mysterious to me, but we have them in Maryland.

https://geodata.md.gov/imap/rest/services/Agriculture/MD_NutrientManagementSetbacksFromWaterways/MapServer/5

"Description: This is a 35 foot buffer area along PDA ditches where fertilizer applications are restricted on croplands."

A PDA is a "Public Ditch Association." Essentially we've already done the environmental damage cutting down the forests and then terraforming the soil for industrial crop harvests. The ditches stretch for miles in some Eastern Shore counties. Taxes pay for the maintenance, as roadside ditches will fill in with dirt and trash otherwise. 

Initially I found maps showing the ditches, and the surrounding service area, as online maps.



And a map of fire coverage.


One good starting point for "ditches" is https://data.imap.maryland.gov/search?collection=dataset&q=ditch. I followed QGis tutorial lessons [https://docs.qgis.org/3.34/en/docs/training_manual/forestry/stands_digitizing.html], digitizing the above into 2 layers, one for the boundaries and one for the lines. Not very good resolution when you zoom in.






I found, after searching with different criteria, similar tax ditch geo-data in the Maryland State site, and the Eastern Shore collective site. So instead of trying to redo the digitizing effort I could skip that part. 

When I downloaded the shape files, they looked good. But when I tried to pull them into a PostgreSQL database, the points flew somewhere else. I am unsure what was missed. I made a second attempt by downloading a KML file instead of a set of shape files. These lines and areas worked, after a fashion. My PostGIS skills are still pretty fresh, like only this year did I run a local database. 


I experimented with arrows for directions, as this helps visualize the land slope, even if miniscule. There is a great example here: https://docs.qgis.org/3.34/en/docs/training_manual/vector_analysis/network_analysis.html. The "loop" above is an oddity from the ditch crossing the highway intersection. Changing the scale and map view sometimes triggers taffy-pull looking connection lines.

This page has more about rules: https://plugins.qgis.org/planet/tag/thematics/.

Next step for me is altering the labels so that roads are distinct from ditches, because there is both a Chicken Bridge Road and a Chicken Bridge drainage area/network. 


Tuesday, February 20, 2024

Well Sample Data and GIS Development

 I looked at water quality tests for a few private wells in Maryland to be able to better use GIS tools. After reviewing public information I took note of key fields such as well depth. In order to show connections among the various aspects I started an entity relationship diagram, going from basic SQL to end up with 8 tables.


In order to keep consistent a few constraints were included, such as:


ALTER TABLE

  well

ADD CONSTRAINT

  well_fkey_property

FOREIGN KEY

  (property_fk)

REFERENCES

  property(name)

;


  I decided to connect wells to a property, and buildings likewise. Normally one well is allowed per property, though while a replacement is being drilled one property might have 2 wells, and if a well runs dry you might say the property has no well.

 The old forms (pre-2000?) have Maryland grid coordinates, with accuracy to 1,000 feet, while the newer forms have a place for latitude/longitude. Water quality tests done by commercial lab contain values of constituents ranging from coliform bacteria to lead and arsenic.

 Sometimes the applications have typos.




Oops, mixed up east and north. I did the same when I first tried to load point data into GIS; once I noticed the drift and switched the values the wellhead came online where it should be.

 After the table definitions and test data load I used the QGis example tutorials to include the bits required for my hand-curated data to be mapped correctly.

See: 
https://docs.qgis.org/3.28/en/docs/training_manual/spatial_databases/simple_feature_model.html
https://docs.qgis.org/3.28/en/docs/training_manual/spatial_databases/geometry.html

> insert into geometry_columns values ('','public','well','the_geometry',2,4326,'POINT');

> update well set the_geometry = 'SRID=4326;POINT(-75.882996 38.938629)' where name = 'CO-15-0020' ;

 






An 8-way join, or is 7-way?

select
  r.value ,                             -- value
  r.uom,                                -- units
  p.name as name_p,                     -- param
  substr(r.name,1,12) as name_r,        -- result
  s.well_fk as well,                    -- well
  s.name as name_s,                     -- sample
  substr(s.sample_place,1,4) as place,  --
  substr(p.descr,1,10) as param_desc,   --
  substr(c.name,1,8) as name_c,         -- corp
  w.depth,                              -- well
  w.pump_installed as pump,             -- well
  substr(o.lot,1,10) as lot,            -- prop
  e.name as person,                     -- people
  s.sample_date
from
  people e,
  building b,
  lab_corp c,
  lab_sample s,
  lab_sample_results r,
  parameter p,
  property o,
  well w
where
  e.building_fk = b.name     -- people in building
and
  b.property_fk = o.name     -- building on property
and
  s.well_fk = w.name         -- sample from well
and
  w.property_fk = o.name     -- well on property
and
  r.parameter_fk = p.name
and
  r.lab_sample_fk = s.name
and
  s.lab_corp_fk = c.name
;


What have I accomplished? Now I have environmental data in a reachable database, which I can edit via command line, or even using LibreOffice forms, and now view and update via QGis.
Next I will figure out the step up from points data (the well heads) to polygons (property lines and building walls). Then, with sample data, set up views with contaminant levels on the map.



Looks a lot like dBASE IV, right?