代做GGR273 Summer 2024 – Lab 3 Health Service Planning调试Python程序

GGR273 Summer 2024 – Lab 3

GGR273 Lab 3: Health Service Planning – PCI Clinic Accessibility in Ontario Health West Region

Due: Aug 12th @ 11:59 pm EST through the Quizzes tab Submit through Lab 3 and answer questions, and submit the final map as a JPEG file.

Objective: Analyze the accessibility of PCI clinics in the Ontario Health West region using QGIS. Perform network analysis, and create a choropleth map showing changes in drive time. Recommend a new PCI facility location.

Context:

PCI (Percutaneous Coronary Intervention) Clinics: PCI clinics are specialized medical facilities equipped to perform. percutaneous coronary interventions, commonly known as angioplasty. This procedure is used to open clogged heart arteries, typically involving the insertion of a balloon-tipped catheter to enlarge the artery, and often the placement of a stent to keep it open. PCI is essential in treating coronary artery disease, preventing heart attacks, and improving overall heart function.CorHealth Ontario is an organization of experts that directs cardiac and stroke services in Ontario. This organization recommends that all residents of Ontario have access to a PCI clinic within a 60-minute drive. This is because from the time that symptoms are apparent, a patient  has a window of 60-minutes to receive care before significant damage is done to their body. The goal of this tutorial is to understand which communities in the Ontario Health West region are outside of this guideline drive time and identify the hospital that does not have a PCI clinic that should be recommended for this investment.

Datasets (download directly from the Quercus modules section as a .zip folder):

1.  Ontario Ministry of Health List of Facilities (provider locations):

https://geohub.lio.gov.on.ca/datasets/d67220ebba164afa9cf38c0f525f2c0a/explore

(original file and reference file)

2.  Population Count by Census Subdivision (CSD):

https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download- telecharger.cfm?Lang=E

3.  CSD Boundary File:https://www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/index2021-eng.cfm?year=21

4.  Ontario Health Regions Boundaries:https://data.ontario.ca/dataset/ontario-s- health-region-geographic-data

Step 1: Importing and Preparing Data

Download the zipped folder of data for Lab 3 from Quercus under Module 6 and 7.

Check the data and ensure that any shapefiles that are in a .zip folder are extracted so they are accessible.

Preparing and Cleaning the Datasets:

1.  Load the files into QGIS

a.  You can click, drag, and drop files into QGIS from the file explorer

directly. To drag and drop shapefiles, select the .shp file. For adding files through the GUI, follow steps b onwards.

b.  You can load data by clicking on the Layer menu drop down and follow the selection path Layer > Add Layer > Add Vector Layer (or “Add

Delimited Text Layer” for the .csv files)

c.  In the Data Source Manager dialog that appears, ensure that the

“Vector” tab is selected (or “ Delimed Text” for the .csv files). For Source Type, choose “ File”

d.  Click the “ …” button next to the file name/datasets field.

e.  In the file browser dialog, navigate to the location of the file and select to add it

f.   For the .csv files: ensure the correct delimiter (commas) is selected under the “ File Format” section.

i.  Verify that the “ First record has field names” checkbox is selected to include the first column as variable headers

ii.  Specify the lat/long fields as X = Longitude and Y = Latitude under the “Geometry Definition” section. If you are loading a table without lat/long you can select “ No Geometry (attribute only table).

2.  Hospital Data (Ministry of Health Service Provider Locations table):

a.  Add a new field called “ PCI_Clinic” as a text field, and for the following hospitals by name (ENGLISH_NA), set them equal to Yes:

The following hospitals have a PCI clinic in Ontario according to the CorHealth website:

Hamilton Health Sciences – General, Health Sciences North - Ramsay Lake Health

Centre, Kingston Health Sciences Centre - Kingston General, London Health Sciences   Centre – Victoria, Peterborough Regional Health Centre, Scarborough Health Network –

Centenary, St. Mary's General Hospital, Sunnybrook Health Sciences Centre - Bayview Campus, Unity Health Toronto - St. Michael's, Southlake Regional Health Centre,

Thunder Bay Regional Health Sciences Centre, Trillium Health Partners- Mississauga, University Health Network - Toronto General, University of Ottawa Heart Institute,

William Osler Health System – Civic, Windsor Regional Hospital - Ouellette Campus, Niagara Health System - St. Catharines

To achieve this:

•   Toggle on editing mode by clicking the pencil icon that is to the farthest left of the Attribute Table ribbon

•   Open the Field Calculator to Create a New Field

•   Set the output name, type and length.

•    In the Expression section, enter the following:

CASE

WHEN "ENGLISH_NA" IN (

'Hamilton Health Sciences - General',

'Health Sciences North - Ramsay Lake Health Centre', 'Kingston Health Sciences Centre - Kingston General',

'London Health Sciences Centre - Victoria', 'Peterborough Regional Health Centre',

'Scarborough Health Network - Centenary', 'St. Mary\'s General Hospital',

'Sunnybrook Health Sciences Centre - Bayview Campus', 'Unity Health Toronto - St. Michael\'s',

'Southlake Regional Health Centre',

'Thunder Bay Regional Health Sciences Centre', 'Trillium Health Partners- Mississauga',

'University Health Network - Toronto General', 'University of Ottawa Heart Institute',

'William Osler Health System - Civic',

'Windsor Regional Hospital - Ouellette Campus',

'Niagara Health System - St. Catharines' ) THEN 'Yes'

ELSE '' END

a.  Save Layer edits and toggle editing off

Step 2: Query the datasets

b.  Filter the service provider locations by service type (SERVICE_TY) and select records labeled as either “ Hospital – Corporation” or “ Hospital – Site). These two different labels identify which hospital is the main address for the corporation of the collection of hospitals for the region that the Ministry of Health uses for funding. Once you filter the records, export it to a new table.

i.         In the attribute table, click on the Select/filter features using Form icon (that looks like a funnel). The exact syntax for each of the two labels is required to filter. You can right click on the cell and copy the cell contents  for both label types. Copy and paste this into the expression box for the SERVICE_TY field and Apply.

ii.          Another way to filter is to click the Advanced Filter (Expression) button in the bottom left corner and enter the following expression. Be sure to double check your syntax by copying the cell content from the labelled cells you want to select:

("SERVICE_TY"='Hospital - Site' OR "SERVICE_TY"='Hospital -    Corporation')

c.  Ensure you have the records selected where SERVICE_TY is hospitals that are labelled as both site and corporation. You can do this by filtering and you should have the option to Export. If not, also select with your cursor by holding down the key to select multiple records that is for your machine (like SHIFT for PC).

d.  Right click on the Layer for the hospitals table and export the selected features as a new shapefile

e.        You should now be working with a file that only has hospitals (labeled as site or   corporation) and an additional column that indicates if the hospital as a PCI clinic. Double check the number of clinics.

f.      Ensure all of your shapefiles and the project are in the projection UTM Zone 17 N

3.  Ontario Health Regions Shapefile

a.  Load this shapefile

b.  Reproject it to UTM Zone 17N

c.  Select the Ontario Health West region and export it so it is a standalone shapefile

d.  We will use this West Region layer to clip the population data so we can focus only on communities that are in the Ontario Health West   region. We want to use hospitals that are for the entire province

because patients can travel to their closest facility, despite Ontario Health boundaries.

4.  Population Data (Population table .csv and CSD shapefile)

The population data include a .csv table that has Census Subdivision (CSD) codes as unique identifiers of each CSD (‘Geography’) as well as the population count from the 2021 Census of Canada. Other data variables are defined in the data dictionary that can be found on the Statistics Canada website. The shapefile of CSD boundaries that matches this file is also included. Join these tables using the CSD unique ID as the primary key then clip the layer to the Ontario Health West region. Ensure all layers are in the   UTM Zone 17N map projection.

Cleaning the population data:

1.  Load the .csv of population data into QGIS

2.  Open the attribute table and review the dataset. Check for the cells that may have values other than numbers. We want to have standard data for the CSD primary key and the population count.

3.  Toggle on editing for this population data table

4.  In the Field Calculator, create a new field called CSD_ID that is the same  data type as the CSD ID field in the shapefile of CSD boundaries that you will join to.

5.  Use the following SQL expression to check for valid CSD IDs and populate the new field with valid CSD IDs. We can see that a valid ID is 7 numbers   long. Anything else will be set to a NULL value:

CASE

WHEN regexp_match("Geography", '^\\d{7}$') THEN "Geography" ELSE NULL

END

6.  Create a new field for the population count data called Pop2021 that is

numeric and only accepts values from the original column (“ Population,

2021) that matches a number format, and the rest of the values are set to  NULL. This will allow the column to be numeric. If any values are copied that arent numeric, the column data type wont be set as such. Ensure you update the fields to match the field names:

CASE

WHEN regexp_match("Population", '^[0-9]+$') THEN to_int("Population") ELSE NULL

END 

7.   Save your edits and toggle off editing

8.   Join the population count data to the CSD boundary shapefile

a.   Right click on the CSD boundary

shapefile. Properties > Joins > Add Join   button (green plus sign in the bottom left)

9.    Join Layer – set to population data .csv file

10. Join Field – set to the CSD ID in the .csv file

11.Target Field – set to the CSDUID in the CSD boundary shapefile. Click OK then OK again

12.Ensure your layer is in the CRS UTM Zone 17N 13.Check the geometry validity of your new layer:

a.  Vector > Geometry Tools > Check Validity

i.   Enter the CSD boundary file, set the output file, and Run to fix geometries

ii.  This means that the layer might be missing spatial

components or it may be corrupted. This tool fixes these issues.

14.   A spatial index is a data structure that allows for queries. The CSD boundary layer doesn’t have this, which can make tools like Clipping

take a long time or even fail. We can create a spatial index:

a.  Right click on the layer > Properties > Source tab

i.   Under the source tab find the Geometry section. There will be a button for Create Spatial Index. Select this.

15.Clip your projected layer of CSD boundaries with population counts to the OH West region using the Clip tool (Vector > Geoprocessing > Clip)

`         Step 3: Create CSD Centroids

16.Now we want to convert our CSD polygons to points called centroids. This  will provide the origin (starting destination) or incident point for our network analysis. We want to calculate the drive distance from each CSD centroid to the nearest PCI centre, this will give us a reasonable approximation of the drive distance for the entire CSD polygon.

17.          This tool can be found in Vector > Geometry Tools > Centroids

18.          Set the Input layer to the polygon layer

19.         Set the output layer to your project folder and click run

20.         Ensure the resulting layer is in the correct CRS

Step 4: Network Analysis

Preparing for Network Analysis by loading a network dataset:

1.  Go to the Geofabrik website to download all features in OpenStreetMap, which will include the road network features. The file you will download is a .pbf file.   This is OpenStreetMap native protocol binary format. It’s a vector format that

provides access to all OSM features.

2.  Navigate to Ontario through the nested structure of geography > North America > Canada > Ontario. Download the file ontario-latest.osm.pbf

3.  Run the following python code to convert the OSM file to a shapefile called roads and load it into your project:

import os

from qgis.core import QgsVectorLayer, QgsProject

# Define the path to the .osm.pbf file and the output directory. Update to your

#path

input_pbf = "C:/data/ontario-latest.osm.pbf" output_dir = "C:/data/output"

output_shapefile = os.path.join(output_dir, "roads.shp")

# Check if output directory exists, create if it doesn't if not os.path.exists(output_dir):

os.makedirs(output_dir)

# Set up the command for ogr2ogr to convert the OSM file to a shapefile

command = f'ogr2ogr -f "ESRI Shapefile" "{output_shapefile}" "{input_pbf}" lines - where "highway IS NOT NULL"'

# Run the command

conversion_result = os.system(command)

# Check if the conversion was successful

if conversion_result == 0 and os.path.exists(output_shapefile):

print("Conversion successful. The shapefile has been created.")

# Load the shapefile into QGIS

layer = QgsVectorLayer(output_shapefile, "Road Network", "ogr")

# Check if the layer is valid if not layer.isValid():

print("Failed to load the layer!") else:

# Add the layer to the QGIS project

QgsProject.instance().addMapLayer(layer)

print("Shapefile loaded successfully and added to the project.")

# Check if necessary fields are present (e.g., length or cost)

field_names = [field.name() for field in layer.fields()] if "length" in field_names or "cost" in field_names:

print("Layer has the necessary fields for network analysis.") else:

print("Layer does not have fields for network analysis")

else:

print("Conversion failed or the shapefile does not exist.")

Check the geometry of the roads file and clip it to the Ontario Health Central and West Regions

1.  Run the check validity tool on the new roads file

2.  Update the spatial index

3.  Create a layer that includes the Ontario Health Central and West Regions,

ensure the geometry is valid and that this layer has a spatial index

4.  Use the clip tool and clip the roads

shapefile to the Ontario Health Central and West Regions.

5.  This layer is now our road network service layer

Step 5: Analyzing Accessibility

1.  Create a network model of CSD centroids to hospitals

a.  Merge your CSD centroids layer and hospital point locations into a single layer, using the Vector > Data Management Tools > Merge   Vector Layers tool.

b.  Create a new field called Name, which will be equal to the CSDUID

field for the points that came from the CSDUID layer, and the

ENGLISH_NA field for the points that came from the hospital layer.

Youll later set this as the unique field that will population the origin and destination fields of the output layer.

c.  Ensure the QNEAT3 plugin is installed

d.  Open QNEAT3 OD Matrix from Points as Lines Tool: Processing

Toolbox > QNEAT3 > OD Matrix from Points as Lines. e.  Configure the tool:

i.   Input point layer: Select your merged points layer containing both CSD and hospital points.

ii.   ID field: Select a unique ID field for the points: likely FID but you can check using the following python code:

layer = QgsProject.instance().mapLayersByName('YourLayerName')[0]

#the following is one line of code

unique_field = next((field.name() for field in layer.fields() if

len(layer.dataProvider().uniqueValues(layer.fields().indexOf(field.name()))) == layer.featureCount()), None)

#the following is one line of

print(f"The field '{unique_field}' has unique values for all features." if unique_field else "No field with unique values was found.")

iii.   Network layer: Select your road network layer

iv.  Travel cost: Choose the appropriate cost field, typically length or distance.

v.  Generated matrix geometry style. matrix geometry follows routes

vi.   Entry Cost Calculation method: Planar

1.  This assumed a flat 2D plane that considers the CRS.  This is fine for the geographic extent of the study area.

vii.  Output layer: Choose the output layer name, which will be lines representing the shortest paths.

viii.  Choose Run

2.  Calculate drive distance information

a.  Consider the list of questions

b.  Use the tools from this course for data handling (Attribute table tools  like Field Calculator, python, and spreadsheet/data tables) to answer the questions at the end of this lab assignment.

c.  You now have the drive distance between all CSD centroids in the Ontario Health West region as well as from all CSD centroids and  hospitals. You also have a field that notes which hospitals are PCI clinics according to CorHealth.

d.  Consider:

i.  Whats the average drive distance (in km) from the CSD centroid to the closest hospital for all hospitals?

ii.  What is the average drive distance from the CSD centroids to

the closest hospital with a PCI clinic? iii.

Step 6: Creating the Choropleth Map

Visualize the drive distance from each CSD centroid in the Ontario Health West region to the closest PCI clinic as a choropleth map. Your final map should have the drive distance layer (in km), location of PCI clinics labelled, location of all hospitals symbolized differently from the PCI clinics, the Ontario Health region boundaries, and   labels for surrounding Great Lakes and the USA border. Verify that all maps elements   are placed and visible: title, legend, scale bar, north arrow, neatline, author name, map authorship date, and data credits. Credits can be included as the organization name followed by a comma, then the year the data was published. Source text can be small and placed to the side of the map as its not the main focus. Your map should be exported as a .jpeg, saved with the course code, lab number and your name. Save your final map as A4 (landscape) and 200 dpi resolution. Open the file to ensure it saved properly before you submit it. Your map will be graded based on the 6 principles of cartography.

Step 7: Answering Questions (20%)

Submit your answers to the following questions in the Quizzes tab by the due date:

1.  (10 points/10%): File submission - submit your map as a .jpeg file. If you do not submit a .jpeg file this section will receive zero points. This file (and all files) will only be accepted over Quercus. 6 points are for each of the 6 principles of cartography. 1 point is for having the required layers listed above and 1 point is for having the correct drive distance data. 2 points are for having the correct symbolization and legend formatting for the drive distance data layer.

2.  (3 points/3%) What is the average drive distance (in km) to the nearest hospital  with a PCI clinic for all CSD centroids in the Ontario Health West region, without adding any new PCI Clinics? Round your answer to the nearest whole number.

3.  (3 points/3%) What would be the average drive distance from the Ontario Health  Region West CSD Centroids to the closest PCI clinic if all hospitals in the Ontario Health West region had a PCI clinic? Round your answer to the nearest whole number.

4.  (1 points/1%) What is the primary key for the CSD layer?

o  Geometry

o  Population, 2021

o  CSD ID

o  Community

5.  (1 points/1%) The length of the new PCI_Clinic field is 100 because it accommodates the range of names of hospitals. True or False.

6.  (1 points/1%) In Field Calculator when setting a new field as equal to a string- type field, for all cells to be equal to the same string, it requires single quotes  around the alphanumeric text. True or False.

7.  (1 points/1%) What CRS does the MMQGIS geocoding service produce a file in when geocoding addresses using the OpenStreetMap reference?

o  EPSG: 32189 – North American Datum 1983 Ontario MNR Lambert

o  EPSG: 4269 – North American Datum 1983

o  EPSG: 3857 – World Geodetic System 1984 Pseudo Mercator

o  EPSG: 4326 – World Geodetic System 1984

o  EPSG: 26917 – North American Datum 1983 – UTM Zone 17N

8.  Submit your QGIS project as a zip file (.zip) using the following instructions (not graded but required to receive full marks. We are using this to ensure you are

completing your own work). Without submitting this, you will receive an initial

point deduction of 5 points and we will follow up to request you to submit this file. This file (and all files) will only be accepted over Quercus.

o  Save all of your QGIS project files and layers.

o  Compile the following into a single .zip file:

1.  QGIS project file as a .qgz file

2.  All shapefiles created during the lab as full sets of sub-files

3.  The exported map file that you submitted in question 1

4.  Atext file/word document that includes the answers to the

questions in this lab with your name and student number on it

Python code in case you get the error in your OD Matrix that your merged file doesn’t have geometry:

from qgis.core import QgsGeometry, QgsPointXY, QgsFeature # Function to validate and update geometry if necessary

def validate_and_update_geometry(feature, default_geometry):

if feature.geometry() is None or not feature.geometry().isGeosValid(): print(f"Updating geometry for feature with ID {feature.id()}")

# Update the geometry to a valid default point (you can customize this as needed) feature.setGeometry(QgsGeometry.fromPointXY(default_geometry))

return True

return False

# Example usage in a loop

default_geometry = QgsPointXY(0, 0) # Default point to use if geometry is invalid # Layer assignment

layer_name = 'Merge_CentroidsHospv3'

layer = QgsProject.instance().mapLayersByName(layer_name)[0] # Iterate through features and update invalid geometries

for feature in layer.getFeatures():

if validate_and_update_geometry(feature, default_geometry):

layer.dataProvider().changeGeometryValues({feature.id(): feature.geometry()})