Plotting points on Australia Map in Python

March 8, 2023

I have the following pandas df

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Define the list of cities and their latitudes/longitudes
cities = ['Sydney', 'Melbourne', 'Brisbane', 'Perth', 'Adelaide', 'Gold Coast', 'Newcastle', 'Canberra', 'Wollongong', 'Sunshine Coast', 'Hobart', 'Alice Springs']
lats = [-33.865143, -37.813629, -27.469771, -31.953512, -34.928499, -28.016666, -32.916668, -35.282001, -34.4278, -26.6500, -42.8821, -23.6970]
longs = [151.209900, 144.963165, 153.025124, 115.857048, 138.600739, 153.399994, 151.750000, 149.128998, 150.8931, 153.0500, 147.3272, 133.8836]

# Generate a random priority for each city between 1 and 5
priority = np.random.randint(1, 6, size=len(cities))

# Create the DataFrame with the city data
data = {'cityname': cities,
        'lats': lats,
        'longs': longs,
        'priority': priority}
df = pd.DataFrame(data)

and I got the shapefile from https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files/STE_2021_AUST_SHP_GDA2020.zip

And the rest of the code are,

# Load the shapefile of Australia
australia = gpd.read_file('STE_2021_AUST_SHP_GDA2020/STE_2021_AUST_GDA2020.shp')

# Define the CRS of the shapefile manually
australia.crs = 'epsg:7844'

# Create a GeoDataFrame from the DataFrame of cities
gdf_cities = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longs, df.lats))

# Set the CRS of the GeoDataFrame to EPSG 7844
# https://epsg.io/7844
gdf_cities.crs = 'epsg:7844'

# Reproject the GeoDataFrame of cities to match the CRS of the shapefile
gdf_cities = gdf_cities.to_crs(australia.crs)

# Perform a spatial join to link the cities to their corresponding polygons in the shapefile
gdf_cities = gpd.sjoin(gdf_cities, australia, predicate='within')

# Set up the plot
sns.set_style("whitegrid")
fig, ax = plt.subplots(figsize=(10, 10))

# Plot the cities colored by priority
sns.scatterplot(ax=ax, data=gdf_cities, x='longs', y='lats', hue='priority', s=100, legend='full')

# Add the shapefile of Australia as a background map
australia.plot(ax=ax, color='lightgrey', edgecolor='white')

# Set the plot title and axis labels
plt.title('Cities in Australia')
plt.xlabel('Longitude')
plt.ylabel('Latitude')

# Show the plot
plt.show()

I could not figure out why the points are not plotted correctly on the Australian map. Can someone review where I am doing it wrong and let me know how to fix? E.g. Alice Springs is in the middle of Australia. Thanks and much appreciated.

>Solution :

Try using the zorder parameter to plot your scatter points on top of the Australian map:

# Plot the cities colored by priority
sns.scatterplot(ax=ax, data=gdf_cities, x='longs', y='lats', hue='priority', s=100, legend='full', zorder=2)

# Add the shapefile of Australia as a background map
australia.plot(ax=ax, color='lightgrey', edgecolor='white', zorder=1)