
import pandas as pd
import numpy as np
import geopandas as gpd # Vector geospatial operations
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import xarray as xr
import imageio # Used for making animated GIFs
from IPython.display import Image
from osgeo import gdal, osr # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
from util import *
plt.rcParams['figure.figsize'] = (20, 20)


# load in pickled geocube 
# loading in a geocube saved as a tif results in a different crs definition (epsg -> wkt) compared to a freshly made geocube, making a comparison of geocubes for consistency fail
import pickle
with open('output/example.pkl', 'rb') as f:
    sample = pickle.load(f)


import sys
sys.executable

'/home/ubuntu/miniconda3/envs/land-use/bin/python'


ls()

Total: 1253.0MB


%%time
REGC = gpd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)

epsg:2193
CPU times: user 2.11 s, sys: 137 ms, total: 2.24 s
Wall time: 15.6 s


land = gpd.read_file('input/lds-nz-coastlines-and-islands-polygons-topo-150k-FGDB.zip!nz-coastlines-and-islands-polygons-topo-150k.gdb')
land['is_land'] = 1


%%time
# TALB['TALB2021_CODE'] = TALB['TALB2021_V1_00'].astype(np.uint32)
cols = ['is_land']
landcube = make_geocube(
    vector_data=land,
    measurements=cols,
    like=sample, # Ensures dimensions match
    fill=0 # non-land will be 0
)
landcube

CPU times: user 5.78 s, sys: 171 ms, total: 5.95 s
Wall time: 5.95 s

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    is_land      (y, x) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[1., 1., 1., ..., 0., 0., 0.],
       [1., 1., 1., ..., 0., 0., 0.],
       [1., 1., 1., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 1., 1., 1.],
       [0., 0., 0., ..., 1., 1., 1.],
       [0., 0., 0., ..., 1., 1., 1.]])


outfile = "output/is_land.tif"
if not os.path.isfile(outfile):
    landcube.astype(np.uint8)["is_land"].rio.to_raster(outfile, dtype=np.uint8)


plt.imshow(landcube.is_land)

<matplotlib.image.AxesImage at 0x7f202658c9a0>


%%time
df = gpd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")

CPU times: user 4min 37s, sys: 20 s, total: 4min 57s
Wall time: 4min 57s


print(df.columns)
print(df.crs)
display(df.sample(5))

Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
       'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
       'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
       'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
       'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
      dtype='object')
epsg:2193


%%time
df = gpd.clip(df, AKL)

CPU times: user 1min 31s, sys: 41.5 ms, total: 1min 32s
Wall time: 1min 31s


urban = [1,2,5]
vegetation = [68,69,71]
water = [0,20,21,22,45,46]


df[df.Class_2018.isin([0])].sample(5)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-b25fb01511d1> in <module>
----> 1 df[df.Class_2018.isin([0])].sample(5)

/data/miniconda3/envs/land_use3/lib/python3.9/site-packages/pandas/core/generic.py in sample(self, n, frac, replace, weights, random_state, axis)
   5348             )
   5349 
-> 5350         locs = rs.choice(axis_length, size=n, replace=replace, p=weights)
   5351         return self.take(locs, axis=axis)
   5352 

mtrand.pyx in numpy.random.mtrand.RandomState.choice()

ValueError: a must be greater than 0 unless no samples are taken


df[df.Name_2018 == 'Mangrove'].sample(5)


for year in years:
    print(year, len(df[df[f'Class_{year}'] == 0]))

1996 3
2001 3
2008 3
2012 3
2018 0


df.Name_2018.value_counts()

Exotic Forest                                3981
Indigenous Forest                            3673
Manuka and/or Kanuka                         2282
Broadleaved Indigenous Hardwoods             1788
Built-up Area (settlement)                   1350
High Producing Exotic Grassland              1326
Mangrove                                     1151
Urban Parkland/Open Space                    1099
Estuarine Open Water                          441
Orchard, Vineyard or Other Perennial Crop     436
Short-rotation Cropland                       362
Lake or Pond                                  326
Herbaceous Saline Vegetation                  303
Low Producing Grassland                       291
Gorse and/or Broom                            287
Forest - Harvested                            266
Sand or Gravel                                252
Deciduous Hardwoods                           201
Surface Mine or Dump                          132
Mixed Exotic Shrubland                        120
Herbaceous Freshwater Vegetation              118
Transport Infrastructure                      107
River                                          15
Flaxland                                        9
Gravel or Rock                                  9
Fernland                                        2
Matagouri or Grey Scrub                         1
Name: Name_2018, dtype: int64


urban = [1,2,5]
vegetation = [68,69,71]
water = [0,20,21,22,45,46]

def simplify_classes(code):
    if code in urban:
        return 1, "Urban"
    elif code in vegetation:
        return 2, "Vegetation"
    elif code in water:
        return 3, "Water"
    else:
        return 4, "Other"

summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    print(year)
    name_year = f"Name_{year}"
    class_year = f"Class_{year}"
    df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
    df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
    display(df[[name_year, class_year]][df[class_year].isin(urban)].value_counts())
    display(df[[name_year, class_year]][df[class_year].isin(vegetation)].value_counts())
    display(df[[name_year, class_year]][df[class_year].isin(water)].value_counts())
    summary.append(df[class_year + "_simplified_name"].value_counts())

1996

Name_1996                   Class_1996
Urban Parkland/Open Space   2             1248
Built-up Area (settlement)  1              489
Transport Infrastructure    5               67
dtype: int64

Name_1996            Class_1996
Exotic Forest        71            3955
Indigenous Forest    69            3737
Deciduous Hardwoods  68             202
dtype: int64

Name_1996                         Class_1996
Estuarine Open Water              22            448
Herbaceous Saline Vegetation      46            313
Lake or Pond                      20            285
Herbaceous Freshwater Vegetation  45            117
River                             21             15
Not land                          0               3
dtype: int64

2001

Name_2001                   Class_2001
Urban Parkland/Open Space   2             1167
Built-up Area (settlement)  1              748
Transport Infrastructure    5               74
dtype: int64

Name_2001            Class_2001
Exotic Forest        71            4495
Indigenous Forest    69            3717
Deciduous Hardwoods  68             206
dtype: int64

Name_2001                         Class_2001
Estuarine Open Water              22            446
Herbaceous Saline Vegetation      46            312
Lake or Pond                      20            305
Herbaceous Freshwater Vegetation  45            116
River                             21             15
Not land                          0               3
dtype: int64

2008

Name_2008                   Class_2008
Built-up Area (settlement)  1             1082
Urban Parkland/Open Space   2             1080
Transport Infrastructure    5              108
dtype: int64


pd.DataFrame(summary).plot.area()

<AxesSubplot:>


%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
    ims = []
    years = [1996, 2001, 2008, 2012, 2018]
    for year in years:
        ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
        ax.set_title(year)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("land_use.gif", ims, fps=1)


with open('land_use.gif','rb') as file:
    display(Image(file.read()))


cols = [f"Class_{year}_simplified_code" for year in years]
cols

['Class_1996_simplified_code',
 'Class_2001_simplified_code',
 'Class_2008_simplified_code',
 'Class_2012_simplified_code',
 'Class_2018_simplified_code']


%%time
geocube = make_geocube(
    vector_data=df,
    output_crs="epsg:2193",
    measurements=cols,
    resolution=(-100, 100),
    fill=0,
)
geocube

CPU times: user 1min 17s, sys: 362 ms, total: 1min 18s
Wall time: 1min 18s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2001_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2008_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2012_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2018_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


# pickle so there is a consistent geocube to compare to later, when producing further data
# (saving and loading a geocube will change crs definition from epsg to wkt, making a comparison of geocubes for consistency fail)
import pickle
with open('output/example.pkl', 'wb') as f:
    pickle.dump(geocube.Class_1996_simplified_code, f, protocol=pickle.HIGHEST_PROTOCOL)

import pickle
with open('output/example.pkl', 'rb') as f:
    sample = pickle.load(f)


# Replace offshore areas (0) with 3 (Water)
geocube = geocube.where(geocube != 0, 3)


geocube.Class_2018_simplified_code.plot()

<matplotlib.collections.QuadMesh at 0x7fb52538fc40>


for year in years:
    print(year)
    outfile = f"output/land_use_{year}.tif"
    if not os.path.isfile(outfile):
        geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize

1996
2001
2008
2012
2018


# check what unique class number / use name pairs there are across all years
years = [1996, 2001, 2008, 2012, 2018]
uses = []
for year in years:
    unique_uses_df = df[[f'Name_{year}', f'Class_{year}']].drop_duplicates()
    uses.extend(unique_uses_df.apply(lambda x: f"{x[f'Class_{year}']:03d}" + ' ' + x[f'Name_{year}'], axis=1).to_list())
# omit duplicates; any use numbers are paired with different use names, or have names spelt differently, they won't be dropped
uses = list(set(uses))
uses.sort(key = lambda x: int(x[:3]))
uses

['000 Not land',
 '001 Built-up Area (settlement)',
 '002 Urban Parkland/Open Space',
 '005 Transport Infrastructure',
 '006 Surface Mine or Dump',
 '010 Sand or Gravel',
 '016 Gravel or Rock',
 '020 Lake or Pond',
 '021 River',
 '022 Estuarine Open Water',
 '030 Short-rotation Cropland',
 '033 Orchard, Vineyard or Other Perennial Crop',
 '040 High Producing Exotic Grassland',
 '041 Low Producing Grassland',
 '045 Herbaceous Freshwater Vegetation',
 '046 Herbaceous Saline Vegetation',
 '047 Flaxland',
 '050 Fernland',
 '051 Gorse and/or Broom',
 '052 Manuka and/or Kanuka',
 '054 Broadleaved Indigenous Hardwoods',
 '056 Mixed Exotic Shrubland',
 '058 Matagouri or Grey Scrub',
 '064 Forest - Harvested',
 '068 Deciduous Hardwoods',
 '069 Indigenous Forest',
 '070 Mangrove',
 '071 Exotic Forest']


# map class names to successive integers (no gaps)
old_use_codes = {u[4:]: int(u[:3]) for i, u in enumerate(uses)}
old_use_codes['Ocean'] = -1
old_use_codes

{'Not land': 0,
 'Built-up Area (settlement)': 1,
 'Urban Parkland/Open Space': 2,
 'Transport Infrastructure': 5,
 'Surface Mine or Dump': 6,
 'Sand or Gravel': 10,
 'Gravel or Rock': 16,
 'Lake or Pond': 20,
 'River': 21,
 'Estuarine Open Water': 22,
 'Short-rotation Cropland': 30,
 'Orchard, Vineyard or Other Perennial Crop': 33,
 'High Producing Exotic Grassland': 40,
 'Low Producing Grassland': 41,
 'Herbaceous Freshwater Vegetation': 45,
 'Herbaceous Saline Vegetation': 46,
 'Flaxland': 47,
 'Fernland': 50,
 'Gorse and/or Broom': 51,
 'Manuka and/or Kanuka': 52,
 'Broadleaved Indigenous Hardwoods': 54,
 'Mixed Exotic Shrubland': 56,
 'Matagouri or Grey Scrub': 58,
 'Forest - Harvested': 64,
 'Deciduous Hardwoods': 68,
 'Indigenous Forest': 69,
 'Mangrove': 70,
 'Exotic Forest': 71,
 'Ocean': -1}


# map class names to successive integers (no gaps)
new_use_codes = {u[4:]: i for i, u in enumerate(uses)}
new_use_codes['Ocean'] = -1
new_use_codes

{'Not land': 0,
 'Built-up Area (settlement)': 1,
 'Urban Parkland/Open Space': 2,
 'Transport Infrastructure': 3,
 'Surface Mine or Dump': 4,
 'Sand or Gravel': 5,
 'Gravel or Rock': 6,
 'Lake or Pond': 7,
 'River': 8,
 'Estuarine Open Water': 9,
 'Short-rotation Cropland': 10,
 'Orchard, Vineyard or Other Perennial Crop': 11,
 'High Producing Exotic Grassland': 12,
 'Low Producing Grassland': 13,
 'Herbaceous Freshwater Vegetation': 14,
 'Herbaceous Saline Vegetation': 15,
 'Flaxland': 16,
 'Fernland': 17,
 'Gorse and/or Broom': 18,
 'Manuka and/or Kanuka': 19,
 'Broadleaved Indigenous Hardwoods': 20,
 'Mixed Exotic Shrubland': 21,
 'Matagouri or Grey Scrub': 22,
 'Forest - Harvested': 23,
 'Deciduous Hardwoods': 24,
 'Indigenous Forest': 25,
 'Mangrove': 26,
 'Exotic Forest': 27,
 'Ocean': -1}


for year in years:
    print(year)
    name_year = f"Name_{year}"
    class_year = f"Class_{year}"
    df[class_year + "_sequential_code"] = df[name_year].apply(lambda c: new_use_codes[c])
df

1996
2001
2008
2012
2018


# either use numbers sequential or the original use codes
# cols = [f"Class_{year}_sequential_code" for year in years]
cols = [f"Class_{year}_sequential_code" for year in years] + [f"Class_{year}" for year in years]
cols

['Class_1996_sequential_code',
 'Class_2001_sequential_code',
 'Class_2008_sequential_code',
 'Class_2012_sequential_code',
 'Class_2018_sequential_code',
 'Class_1996',
 'Class_2001',
 'Class_2008',
 'Class_2012',
 'Class_2018']


%%time
geocube = make_geocube(
    vector_data=df,
    measurements=cols,
    fill=-1, # ocean
    like=sample
)
geocube

CPU times: user 2min 1s, sys: 312 ms, total: 2min 1s
Wall time: 2min 1s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_sequential_code  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2001_sequential_code  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2008_sequential_code  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2012_sequential_code  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2018_sequential_code  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_1996                  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2001                  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2008                  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2012                  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
    Class_2018                  (y, x) float64 -1.0 -1.0 -1.0 ... -1.0 -1.0 -1.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       ...,
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.]])

array([[-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       ...,
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.]])

array([[-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       ...,
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.]])

array([[-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       ...,
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.]])

array([[-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       ...,
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.],
       [-1., -1., -1., ..., -1., -1., -1.]])


# look carefully to see 'Not land' areas - some kind of dock infrastructure, perhaps
# appear to be safe to turn into ocean
year = 2008
geoarray = np.array(geocube[f'Class_{year}'])
geoarray[np.where(geoarray > 0.01)] = 1
plt.figure(figsize=(30,30))
plt.matshow(geoarray[600:1050, 380:630], fignum=1)
plt.show()


for year in years:
    print(year)
    outfile_sequential = f"output/land_use_detailed_sequential_{year}.tif"
    outfile = f"output/land_use_detailed_{year}.tif"
    if not os.path.isfile(outfile):
        # detailed land uses mapped to sequential integers (-1 for ocean)
        geocube[f"Class_{year}_sequential_code"].rio.to_raster(outfile_sequential, dtype=np.byte) # Use np.byte for smaller output filesize
        # detailed land uses with original integers (-1 for ocean)
        geocube[f"Class_{year}"].rio.to_raster(outfile, dtype=np.byte)

1996
2001
2008
2012
2018


%%time
TALB = gpd.read_file("input/statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip!territorial-authority-local-board-2021-clipped-generalised.gdb")
TALB = gpd.clip(TALB, AKL)
TALB.head()

CPU times: user 767 ms, sys: 118 ms, total: 885 ms
Wall time: 877 ms


# From https://www.stats.govt.nz/information-releases/subnational-population-projections-2018base2048#map
pop = pd.concat(pd.read_excel(
    "input/subnational-population-projections-2018base-2048.xlsx",
    sheet_name=["Table 5", "Table 6"],
    skiprows=5,
    usecols="A,B,G",
    names=["area", "year", "population"],
    nrows=231,
    engine='openpyxl'
)).fillna(method='ffill').reset_index()
pop.replace("Maungakiekie-Tamaki local board area", "Maungakiekie-Tāmaki local board area", inplace=True)
pop.year = pop.year.astype(str)
pop


pop[pop.area.str.contains("local")].pivot(index='year', columns='area', values='population').plot()

<AxesSubplot:xlabel='year'>


pop = pop.pivot(index='area', columns='year', values='population')
pop


cols = pop.columns.tolist()
cols

['1996',
 '2001',
 '2006',
 '2013',
 '2018',
 '2023',
 '2028',
 '2033',
 '2038',
 '2043',
 '2048']


pop = pd.merge(TALB, pop, left_on=TALB.TALB2021_V1_00_NAME.str.lower().str.strip(), right_on=pop.index.str.lower().str.strip(), how="left")
pop[["TALB2021_V1_00_NAME"] + cols]


d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()

(17800, 84495.2479338843, 190900)


%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
    ims = []
    for col in cols:
        ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
        ax.set_title(col)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("pop.gif", ims, fps=1)


with open('pop.gif','rb') as file:
    display(Image(file.read()))


%%time
popcube = make_geocube(
    vector_data=pop,
    measurements=cols,
    like=geocube, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
popcube

CPU times: user 1.66 s, sys: 9.95 ms, total: 1.67 s
Wall time: 1.67 s

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    1996         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2001         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2006         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2013         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2018         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2023         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2028         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2033         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2038         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2043         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2048         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


for col in cols:
    outfile = f"output/pop_{col}.tif"
    if not os.path.isfile(outfile):
        # byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
        popcube[col].rio.to_raster(outfile, dtype=np.uint32)


%%time
# create a 'talbcube' also, identifying which areas belong to which talbs
# this will enable easy segmentation into areas at which population is defined, for operations like calculating density
TALB['TALB2021_CODE'] = TALB['TALB2021_V1_00'].astype(np.uint32)
cols = ['TALB2021_CODE']
talbcube = make_geocube(
    vector_data=TALB,
    measurements=cols,
    like=sample, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
talbcube

CPU times: user 3.2 s, sys: 139 ms, total: 3.34 s
Wall time: 3.34 s

<xarray.Dataset>
Dimensions:        (x: 1001, y: 1320)
Coordinates:
  * y              (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x              (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref    int64 0
Data variables:
    TALB2021_CODE  (y, x) float64 300.0 300.0 300.0 ... 1.3e+03 1.3e+03 1.3e+03

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[ 300.,  300.,  300., ...,    0.,    0.,    0.],
       [ 300.,  300.,  300., ...,    0.,    0.,    0.],
       [ 300.,  300.,  300., ...,    0.,    0.,    0.],
       ...,
       [   0.,    0.,    0., ..., 1300., 1300., 1300.],
       [   0.,    0.,    0., ..., 1300., 1300., 1300.],
       [   0.,    0.,    0., ..., 1300., 1300., 1300.]])


for col in cols:
    outfile = f"output/talb_{col}.tif"
    if not os.path.isfile(outfile):
        # byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
        # first using astype seems necessary to avoid scary errors
        talbcube.astype(np.uint32)[col].rio.to_raster(outfile, dtype=np.uint32)


%%time
roads = gpd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")

CPU times: user 41 s, sys: 1.3 s, total: 42.3 s
Wall time: 42.3 s


%%time
akl_roads = gpd.clip(roads, AKL)

CPU times: user 43.4 s, sys: 968 µs, total: 43.4 s
Wall time: 43.4 s


# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway


mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()


mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()


ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)


%%time
mway_cube = make_geocube(
    vector_data=mway,
    measurements=["lane_count"],
    like=sample, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube

CPU times: user 672 ms, sys: 20 ms, total: 692 ms
Wall time: 691 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    lane_count   (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
    mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)


src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


%%time
skytower = [-36.84838748948485, 174.7621736911587]
skytower = gpd.points_from_xy(x=[skytower[1]], y=[skytower[0]])
cbd = gpd.GeoDataFrame([{"name": "Skytower", "value": 1}], geometry=skytower, crs="EPSG:4326").to_crs(epsg=2193)
display(cbd)
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
    vector_data=cbd,
    like=sample, # Ensures dimensions match
    fill=0
)
display(cbd_cube)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
    cbd_cube.value.rio.to_raster(outfile, dtype=np.byte)
    
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

(1320, 1001)
CPU times: user 880 ms, sys: 60.1 ms, total: 940 ms
Wall time: 935 ms

Text(0.5, 1.0, 'Distance (m)')


bounds = AKL.total_bounds.tolist()
bounds

[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]


zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles

['EJ.tif',
 'DM.tif',
 'EL.tif',
 'DL.tif',
 'DJ.tif',
 'FK.tif',
 'DK.tif',
 'EK.tif',
 'FL.tif']


tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
tile_datasets = [tile for tile in tile_datasets if not rasterio.coords.disjoint_bounds(tile.bounds, bounds)]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100))


DEM[DEM == -32767] = np.nan
print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(DEM, cmap='terrain')

-30.0 78.13751 697.73083 (1, 1320, 1018)

<AxesSubplot:>


cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape

(1, 1320, 1001)


warped_DEM = np.full_like(cbd_dist.read(), np.nan, dtype=np.float32)

rasterio.warp.reproject(DEM,
                        warped_DEM,
                        src_transform=transform,
                        src_nodata=np.nan,
                        dst_nodata=np.nan,
                        src_crs=tile_datasets[0].crs,
                        dst_crs=cbd_dist.crs,
                        dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)

(1, 1320, 1018) (1, 1320, 1001)


meta = tile_datasets[0].meta
print(meta)
meta.update({
    "nodata": np.nan,
    "height": cbd_dist.height,
    "width": cbd_dist.width,
    "transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
    with rasterio.open(outfile, "w", **meta) as dest:
       dest.write(warped_DEM)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
       0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
       0.0, -100.0, 6002400.0)}


gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
slope[slope == -9999] = np.nan
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(slope, cmap='terrain')

0.0 5.5326476 52.67306 (1, 1320, 1001)

<AxesSubplot:>


consents = pd.read_csv("input/new-dwellings-consented-by-statistical-area-2-april-2021-csv.zip")
consents


consents["year"] = pd.to_datetime(consents.month).dt.year.astype(str)
consents = consents.groupby(["year", "SA2_code"]).total_dwelling_units.sum().reset_index()
consents = consents.pivot(index='SA2_code', columns='year', values='total_dwelling_units')


sa2 = gpd.read_file("input/statsnzstatistical-area-2-2021-generalised-FGDB.zip!statistical-area-2-2021-generalised.gdb")


%%time
sa2 = gpd.clip(sa2, AKL)

CPU times: user 2.57 s, sys: 0 ns, total: 2.57 s
Wall time: 2.57 s


sa2.SA22021_V1_00 = sa2.SA22021_V1_00.astype(int)


consents = sa2.merge(consents, left_on="SA22021_V1_00", right_on="SA2_code")


consents


consents.describe()


%%time
years = ["2018", "2019", "2020"]
consentcube = make_geocube(
    vector_data=consents,
    measurements=years,
    like=sample, # Ensures dimensions match
    fill=0 # NaNs will be 0
)
consentcube

CPU times: user 969 ms, sys: 21.9 ms, total: 991 ms
Wall time: 989 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    2018         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2019         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2020         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


consentcube["2018"].plot()

<matplotlib.collections.QuadMesh at 0x7f0ea9a32e80>


for year in years:
    outfile = f"output/consents_{year}.tif"
    if not os.path.isfile(outfile):
        # byte max value is 255, and we have larger values than that here. uint16 max value is 65535
        consentcube[year].rio.to_raster(outfile, dtype=np.uint16)


bus = gpd.read_file('input/BusService.geojson')

bus


bus['value'] = 1
bus = bus[['geometry', 'value']]
bus


%%time
buscube = make_geocube(
    vector_data=bus,
    like=sample, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
display(buscube)
outfile = "output/bus.tif"
if not os.path.isfile(outfile):
    buscube.value.rio.to_raster(outfile, dtype=np.byte)

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

CPU times: user 1.58 s, sys: 17.6 ms, total: 1.59 s
Wall time: 1.59 s


f = "output/bus.tif"
dst_ds = gdal.Open(f)
img_array = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
plt.imshow(img_array)

<matplotlib.image.AxesImage at 0x7fb4e9979760>


src_ds = gdal.Open("output/bus.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/bus_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from bus stop in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
bus_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(bus_dist.shape)
plt.imshow(bus_dist)
plt.title("Distance from Bus Stops in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


train = gpd.read_file('input/TrainService.geojson')


train


train['value'] = 1
train = train[['geometry', 'value']]
train


%%time
traincube = make_geocube(
    vector_data=train,
    like=sample, # Ensures dimensions match
    fill=0,
)
display(traincube)
outfile = "output/train.tif"
if not os.path.isfile(outfile):
    traincube.value.rio.to_raster(outfile, dtype=np.byte)

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

CPU times: user 236 ms, sys: 1.05 ms, total: 237 ms
Wall time: 234 ms


f = "output/train.tif"
dst_ds = gdal.Open(f)
img_array = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
plt.imshow(img_array)

<matplotlib.image.AxesImage at 0x7f0ea8549e20>


src_ds = gdal.Open("output/train.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/train_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from bus stop in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
bus_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(bus_dist.shape)
plt.imshow(bus_dist)
plt.title("Distance from train stations in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


carparks = gpd.read_file('input/CarParking.geojson')


carparks['value'] = 1
carparks = carparks[carparks.DATATYPE.isin(['Open Air Parking', 'Covered Parking'])][['geometry', 'value']]
carparks


%%time
parkcube = make_geocube(
    vector_data=carparks,
    like=sample, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
display(parkcube)
outfile = "output/park.tif"
if not os.path.isfile(outfile):
    parkcube.value.rio.to_raster(outfile, dtype=np.byte)

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

CPU times: user 527 ms, sys: 19.9 ms, total: 547 ms
Wall time: 550 ms


f = "output/park.tif"
dst_ds = gdal.Open(f)
img_array = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
plt.imshow(img_array)

<matplotlib.image.AxesImage at 0x7fb4dbc60670>


src_ds = gdal.Open("output/park.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/park_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from bus stop in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
park_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(park_dist.shape)
plt.imshow(park_dist)
plt.title("Distance from Parks in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


schools = pd.read_csv('input/schools-govt.csv', skiprows=15)
schools = gpd.GeoDataFrame(schools, geometry=gpd.points_from_xy(schools.Longitude, schools.Latitude)).set_crs(4326).to_crs(2193)


ax = schools.sample(100).to_crs(3857).plot()
ctx.add_basemap(ax)


schools['value'] = 1
schools = schools[['geometry', 'value']]
schools


%%time
schoolcube = make_geocube(
    vector_data=schools,
    like=sample, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
display(schoolcube)
outfile = "output/schools.tif"
if not os.path.isfile(outfile):
    schoolcube.value.rio.to_raster(outfile, dtype=np.byte)

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

CPU times: user 599 ms, sys: 20 ms, total: 619 ms
Wall time: 617 ms


f = "output/schools.tif"
dst_ds = gdal.Open(f)
img_array = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
plt.imshow(img_array)

<matplotlib.image.AxesImage at 0x7ff1edd798e0>


src_ds = gdal.Open("output/schools.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/schools_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from bus stop in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
schools_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(schools.shape)
plt.imshow(schools_dist)
plt.title("Distance from Schools in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(559, 2)

Text(0.5, 1.0, 'Distance (m)')


hospitals = gpd.read_file("input/NZ_Public_Hospitals.geojson")
hospitals


hospitals['value'] = 1
hospitals = hospitals[['geometry', 'value']]
hospitals


%%time
hospitalcube = make_geocube(
    vector_data=hospitals,
    like=sample, # Ensures dimensions match
    fill=0,
)
display(hospitalcube)
outfile = "output/hospitals.tif"
if not os.path.isfile(outfile):
    hospitalcube.value.rio.to_raster(outfile, dtype=np.byte)

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

CPU times: user 240 ms, sys: 2.17 ms, total: 242 ms
Wall time: 239 ms


f = "output/hospitals.tif"
dst_ds = gdal.Open(f)
img_array = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
plt.imshow(img_array)

<matplotlib.image.AxesImage at 0x7f0ea94abf70>


src_ds = gdal.Open("output/hospitals.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/hospitals_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from bus stop in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
bus_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(bus_dist.shape)
plt.imshow(bus_dist)
plt.title("Distance from Hospitals in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


check_output()

{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}


[(f, rasterio.open(f"output/{f}").crs) for f in os.listdir("output") if f.endswith(".tif")]

[('land_use_2001.tif', CRS.from_epsg(2193)),
 ('cbd_dist.tif', CRS.from_epsg(2193)),
 ('dem.tif', CRS.from_epsg(2193)),
 ('hospitals_dist.tif', CRS.from_epsg(2193)),
 ('pop_2023.tif', CRS.from_epsg(2193)),
 ('pop_1996.tif', CRS.from_epsg(2193)),
 ('land_use_2012.tif', CRS.from_epsg(2193)),
 ('land_use_2008.tif', CRS.from_epsg(2193)),
 ('schools_dist.tif', CRS.from_epsg(2193)),
 ('pop_2028.tif', CRS.from_epsg(2193)),
 ('pop_2013.tif', CRS.from_epsg(2193)),
 ('park.tif', CRS.from_epsg(2193)),
 ('pop_2018.tif', CRS.from_epsg(2193)),
 ('pop_2006.tif', CRS.from_epsg(2193)),
 ('bus_dist.tif', CRS.from_epsg(2193)),
 ('pop_2001.tif', CRS.from_epsg(2193)),
 ('pop_2038.tif', CRS.from_epsg(2193)),
 ('land_use_1996.tif', CRS.from_epsg(2193)),
 ('bus.tif', CRS.from_epsg(2193)),
 ('pop_2048.tif', CRS.from_epsg(2193)),
 ('pop_2043.tif', CRS.from_epsg(2193)),
 ('cbd.tif', CRS.from_epsg(2193)),
 ('park_dist.tif', CRS.from_epsg(2193)),
 ('consents_2019.tif', CRS.from_epsg(2193)),
 ('hospitals.tif', CRS.from_epsg(2193)),
 ('slope.tif', CRS.from_epsg(2193)),
 ('mway.tif', CRS.from_epsg(2193)),
 ('consents_2018.tif', CRS.from_epsg(2193)),
 ('land_use_2018.tif', CRS.from_epsg(2193)),
 ('mway_dist.tif', CRS.from_epsg(2193)),
 ('schools.tif', CRS.from_epsg(2193)),
 ('consents_2020.tif', CRS.from_epsg(2193)),
 ('pop_2033.tif', CRS.from_epsg(2193))]


ls("output")

Total: 115.0MB

	name	filesize (MB)	last modified
0	.ipynb_checkpoints	0.00	2021-07-01 01:11:29.638198
1	2018-census-main-means-of-travel-to-work-by-statistical-area.csv	6.93	2021-07-06 04:45:38.570015
2	BusRoutes.geojson	13.87	2021-07-06 04:45:39.202016
3	BusService.geojson	1.84	2021-06-17 04:22:11.438204
4	CarParking.geojson	0.08	2021-06-17 05:20:35.495894
5	Crash_Analysis_System_(CAS)_data.zip	47.04	2021-07-06 04:45:39.578016
6	NZ_Public_Hospitals.geojson	0.10	2021-06-29 04:50:34.459662
7	ParkAndRide.geojson	0.01	2021-06-17 05:19:24.551777
8	SA2+IMD+consents+ethnicity.topojson	4.82	2021-06-17 21:34:25.436450
9	TrainService.geojson	0.01	2021-06-17 04:21:04.250090
10	auckland-transport-bus-performance-report-for-01-january-2019-through-30-april-2021.xlsx	3.12	2021-07-06 04:45:38.018014
11	lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip	662.77	2021-06-16 00:07:17.674683
12	lds-nz-road-centrelines-topo-150k-FGDB.zip	35.50	2021-06-16 00:06:56.854655
13	lds-protected-areas-FGDB.zip	44.29	2021-07-21 01:27:04.095702
14	lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip	363.55	2021-06-16 00:07:09.582672
15	new-dwellings-consented-by-statistical-area-2-april-2021-csv.zip	7.60	2021-06-29 04:50:34.587662
16	place_of_worship-OSM.geojson	1.73	2021-07-06 04:45:38.358015
17	schools-OSM.geojson	3.79	2021-06-17 04:52:31.013133
18	schools-govt.csv	0.27	2021-06-17 05:30:21.956860
19	statistical-area-2-2018-generalised.gdb	0.00	2021-06-21 02:27:39.602486
20	statsnzarea-unit-2013-FGDB.zip	13.79	2021-06-18 05:04:42.313101
21	statsnzregional-council-2021-clipped-generalised-FGDB.zip	4.97	2021-06-16 00:06:57.562656
22	statsnzstatistical-area-2-2018-generalised-FGDB.zip	15.47	2021-06-21 02:28:09.406536
23	statsnzstatistical-area-2-2021-generalised-FGDB.zip	14.56	2021-06-29 04:50:34.547662
24	statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip	6.05	2021-06-16 00:06:56.990655
25	subnational-population-projections-2018base-2048.xlsx	0.14	2021-06-16 00:06:58.906658
26	supermarkets-OSM.geojson	0.26	2021-07-06 04:45:38.682015

	TALB2021_V1_00	TALB2021_V1_00_NAME	TALB2021_V1_00_NAME_ASCII	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	geometry
0	00100	Far North District	Far North District	6686.751655	6699.099475	2.167299e+06	MULTIPOLYGON (((1643927.274 6083897.824, 16439...
1	00200	Whangarei District	Whangarei District	2712.118612	2712.118612	9.223848e+05	MULTIPOLYGON (((1752996.594 6015585.829, 17530...
2	00300	Kaipara District	Kaipara District	3108.960347	3108.960347	9.341178e+05	MULTIPOLYGON (((1714357.617 5983211.997, 17143...
3	01100	Thames-Coromandel District	Thames-Coromandel District	2207.766906	2207.766906	8.204000e+05	MULTIPOLYGON (((1857357.606 5876615.882, 18573...
4	01200	Hauraki District	Hauraki District	1270.133527	1270.133527	3.112825e+05	MULTIPOLYGON (((1824948.875 5883405.558, 18249...

year	1996	2001	2006	2013	2018	2023	2028	2033	2038	2043	2048
area
Albert-Eden local board area	86100	90000	96200	100000	103700	106900	111200	116500	121600	126300	130600
Auckland	1115800	1218300	1373000	1493200	1654800	1778700	1891800	2001800	2107000	2207800	2302900
Devonport-Takapuna local board area	51600	52300	55700	58500	60500	62700	65700	68800	71800	74800	77600
Far North district	54500	56400	57500	60600	67900	72400	75000	77300	79100	80600	81700
Franklin local board area	48400	52900	60600	68300	77700	85700	96900	109400	122000	134500	146900
Gisborne district	47200	45500	45900	47000	49500	51900	53100	54000	54600	55000	55200
Great Barrier local board area	1230	1100	940	950	960	1030	1080	1120	1150	1180	1200
Hamilton city	113500	121200	134800	150200	168600	183000	194400	205400	216000	226500	236600
Hauraki district	18550	18000	18300	18600	20700	21800	22200	22400	22300	22100	21800
Henderson-Massey local board area	82200	90900	103700	113500	124600	134100	142500	150500	158200	165700	172900
Hibiscus and Bays local board area	68000	75100	85200	94000	108500	117700	123500	126500	128900	130900	132400
Howick local board area	80900	97300	119100	135000	149400	160500	167800	174400	180400	185900	190900
Kaipara district	17800	17950	18550	20500	23700	25900	27200	28300	29100	29800	30300
Kaipātiki local board area\n	74700	78200	83600	87000	92900	97300	99400	101100	102500	103600	104700
Kawerau district	8120	7290	7150	6650	7460	7910	8000	8020	7970	7860	7720
Manurewa local board area	60100	69100	81600	87000	100900	109600	114000	118000	121700	124900	127600
Matamata-Piako district	30300	30300	31200	32900	35300	37000	37900	38700	39200	39500	39600
Maungakiekie-Tāmaki local board area	61700	66400	70400	73700	80500	86600	93800	101300	108800	116100	123000
Māngere-Ōtāhuhu local board area\n	61400	64600	72500	75300	82700	88700	92500	96300	99800	103100	105900
Papakura local board area	38000	38900	43100	48200	61100	70000	77000	81700	86100	90500	95000
Puketāpapa local board area\n	46100	49800	53900	56300	60900	64600	68800	73000	77100	81000	84600
Rodney local board area	39000	44100	51000	57300	69100	77600	87600	99600	111700	123800	135800
Rotorua district	66600	66900	68100	68400	74800	78900	80700	82200	83400	84200	84800
South Waikato district	25800	24200	23200	23200	24800	25800	26400	26800	27000	27100	27100
Taupō district\n	31600	32500	33400	34800	38600	40900	42000	42800	43300	43600	43800
Tauranga city	79900	93600	107000	119900	142100	156900	166300	175000	183300	191400	199100
Thames-Coromandel district	25400	25800	26700	27300	30700	32400	33100	33500	33500	33300	32800
Upper Harbour local board area	24300	33700	45100	56800	66800	75900	86000	95900	105800	115400	124900
Waiheke local board area	6680	7560	8150	8630	9360	9830	10250	10650	10900	11100	11200
Waikato district	52000	53700	59500	66500	78200	86100	92800	99300	105700	111900	117700
Waipa district	38400	40000	43700	48700	55000	59300	62100	64600	66900	68900	70700
Waitematā local board area\n	44900	52500	66600	81300	88500	94700	102800	111900	120700	129200	137400
Waitomo district	10000	9780	9680	9340	9630	9740	9730	9670	9540	9330	9070
Waitākere Ranges local board area\n	42100	44200	47500	50700	54200	56900	58300	59500	60300	60900	61100
Western Bay of Plenty district	35600	38900	42900	45400	53300	58100	60900	63300	65200	66700	68000
Whakatāne district\n	34200	34100	34500	34200	37100	38800	39300	39500	39500	39300	38900
Whangārei district\n	68400	70000	76500	83700	94100	101000	105500	109500	113100	116400	119300
Whau local board area	59700	66200	73100	76700	84100	89500	95000	100500	105900	111000	115900
Ōpōtiki district\n	9620	9490	9200	8780	9670	10250	10350	10400	10300	10150	9910
Ōrākei local board area\n	70900	73000	78400	83700	87700	91300	96200	101100	105800	110300	114500
Ōtara-Papatoetoe local board area\n	67800	70400	76600	80300	90500	97700	101500	104000	106000	107500	108600
Ōtorohanga district\n	9960	9590	9310	9590	10500	11050	11300	11550	11750	11900	12000

	TALB2021_V1_00_NAME	1996	2001	2006	2013	2018	2023	2028	2033	2038	2043	2048
0	Kaipara District	17800	17950	18550	20500	23700	25900	27200	28300	29100	29800	30300
1	Hauraki District	18550	18000	18300	18600	20700	21800	22200	22400	22300	22100	21800
2	Waikato District	52000	53700	59500	66500	78200	86100	92800	99300	105700	111900	117700
3	Rodney Local Board Area	39000	44100	51000	57300	69100	77600	87600	99600	111700	123800	135800
4	Hibiscus and Bays Local Board Area	68000	75100	85200	94000	108500	117700	123500	126500	128900	130900	132400
5	Upper Harbour Local Board Area	24300	33700	45100	56800	66800	75900	86000	95900	105800	115400	124900
6	Kaipātiki Local Board Area	74700	78200	83600	87000	92900	97300	99400	101100	102500	103600	104700
7	Devonport-Takapuna Local Board Area	51600	52300	55700	58500	60500	62700	65700	68800	71800	74800	77600
8	Henderson-Massey Local Board Area	82200	90900	103700	113500	124600	134100	142500	150500	158200	165700	172900
9	Waitākere Ranges Local Board Area	42100	44200	47500	50700	54200	56900	58300	59500	60300	60900	61100
10	Waitematā Local Board Area	44900	52500	66600	81300	88500	94700	102800	111900	120700	129200	137400
11	Whau Local Board Area	59700	66200	73100	76700	84100	89500	95000	100500	105900	111000	115900
12	Albert-Eden Local Board Area	86100	90000	96200	100000	103700	106900	111200	116500	121600	126300	130600
13	Puketāpapa Local Board Area	46100	49800	53900	56300	60900	64600	68800	73000	77100	81000	84600
14	Ōrākei Local Board Area	70900	73000	78400	83700	87700	91300	96200	101100	105800	110300	114500
15	Maungakiekie-Tāmaki Local Board Area	61700	66400	70400	73700	80500	86600	93800	101300	108800	116100	123000
16	Howick Local Board Area	80900	97300	119100	135000	149400	160500	167800	174400	180400	185900	190900
17	Māngere-Ōtāhuhu Local Board Area	61400	64600	72500	75300	82700	88700	92500	96300	99800	103100	105900
18	Ōtara-Papatoetoe Local Board Area	67800	70400	76600	80300	90500	97700	101500	104000	106000	107500	108600
19	Manurewa Local Board Area	60100	69100	81600	87000	100900	109600	114000	118000	121700	124900	127600
20	Papakura Local Board Area	38000	38900	43100	48200	61100	70000	77000	81700	86100	90500	95000
21	Franklin Local Board Area	48400	52900	60600	68300	77700	85700	96900	109400	122000	134500	146900

	t50_fid	name_ascii	macronated	name	hway_num	rna_sufi	lane_count	way_count	status	surface	geometry
512	100120610	KAIPARA COAST HIGHWAY	N	KAIPARA COAST HIGHWAY	16	3007739	2	None	None	sealed	LINESTRING (1732000.000 5944172.070, 1732048.5...
2933	3198057	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748581.508 5968975.145, 1748558.4...
2934	3198059	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748171.047 5971284.152, 1748129.9...
3320	3200754	PAERATA ROAD	N	PAERATA ROAD	22	3000260	2	None	None	sealed	LINESTRING (1767236.112 5888088.508, 1767244.3...
3324	3200792	UPPER HARBOUR MOTORWAY	N	UPPER HARBOUR MOTORWAY	18	3047073	4	None	None	sealed	LINESTRING (1747954.314 5927269.837, 1747970.0...
...	...	...	...	...	...	...	...	...	...	...	...
138240	100048291	AUCKLAND-WAIWERA MOTORWAY	N	AUCKLAND-WAIWERA MOTORWAY	1	3067966	7	None	None	sealed	LINESTRING (1755881.018 5922863.734, 1755886.4...
138301	100048432	AUCKLAND-HAMILTON MOTORWAY	N	AUCKLAND-HAMILTON MOTORWAY	1	3017109	1	None	None	sealed	LINESTRING (1765115.647 5909916.697, 1765092.7...
138337	100048532	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	4	None	None	sealed	LINESTRING (1748892.089 5949596.727, 1748892.0...
138369	100048589	PORT ALBERT ROAD	N	PORT ALBERT ROAD	16	3013274	2	None	None	sealed	LINESTRING (1734173.019 5980575.187, 1734175.6...
138680	100118365	SOUTH-WESTERN MOTORWAY	N	SOUTH-WESTERN MOTORWAY	20	3018532	4	None	None	sealed	LINESTRING (1760066.252 5908184.133, 1760043.9...

Land use¶

Land Uses Detailed¶

Population density¶

Distance from major roads¶

Distance from CBD¶

Slope¶

Building Consents¶

bus stops¶

Train stations¶

Car Parks¶

Schools¶

Hospitals¶

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
492762	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	71	71	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb2000518016	MULTIPOLYGON (((1510972.655 5181146.580, 15109...
49020	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000186502	MULTIPOLYGON (((1367760.677 5192171.249, 13677...
397826	Low Producing Grassland	Low Producing Grassland	Low Producing Grassland	Low Producing Grassland	Low Producing Grassland	41	41	41	41	41	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000074698	MULTIPOLYGON (((1584254.121 5419266.963, 15841...
446072	High Producing Exotic Grassland	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	40	71	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb1000199223	MULTIPOLYGON (((1808177.512 5570895.577, 18081...
141315	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	54	54	54	54	54	...	no	yes	yes	yes	yes	yes	Landcare Research	2015-06-30T00:00:00	lcdb2000206411	MULTIPOLYGON (((1187228.828 4936074.249, 11872...

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
509195	Mangrove	Mangrove	Mangrove	Mangrove	Mangrove	70	70	70	70	70	...	yes	yes	yes	yes	yes	yes	Landcare Research	2004-06-30T00:00:00	lcdb1000183153	POLYGON ((1759645.787 5904824.637, 1759583.342...
355	Mangrove	Mangrove	Mangrove	Mangrove	Mangrove	70	70	70	70	70	...	yes	no	no	no	no	no	Terralink	2004-06-30T00:00:00	lcdb1000182039	POLYGON ((1747738.916 5923907.316, 1747718.378...
368494	Mangrove	Mangrove	Mangrove	Mangrove	Mangrove	70	70	70	70	70	...	yes	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb1000511804	POLYGON ((1769033.367 5912506.704, 1769057.533...
174	Mangrove	Mangrove	Mangrove	Mangrove	Mangrove	70	70	70	70	70	...	yes	no	no	no	no	no	Terralink	2004-06-30T00:00:00	lcdb1000181804	POLYGON ((1750115.743 5883003.616, 1750132.151...
406436	Mangrove	Mangrove	Mangrove	Mangrove	Mangrove	70	70	70	70	70	...	yes	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000183668	POLYGON ((1718738.141 5978578.636, 1718712.184...

	month	SA2_code	SA2_name	territorial_authority	total_dwelling_units	houses	apartments	retirement_village_units	townhouses_flats_units_other
0	1990-04-01	100100	North Cape	Far North District	2	2	0	0	0
1	1990-04-01	100200	Rangaunu Harbour	Far North District	1	1	0	0	0
2	1990-04-01	100300	Inlets Far North District	Far North District	0	0	0	0	0
3	1990-04-01	100400	Karikari Peninsula	Far North District	2	2	0	0	0
4	1990-04-01	100500	Tangonge	Far North District	1	1	0	0	0
...	...	...	...	...	...	...	...	...	...
842602	2021-04-01	400012	Oceanic Snares Islands	Area Outside Territorial Authority	0	0	0	0	0
842603	2021-04-01	400013	Snares Islands	Area Outside Territorial Authority	0	0	0	0	0
842604	2021-04-01	400014	Oceanic Antipodes Islands	Area Outside Territorial Authority	0	0	0	0	0
842605	2021-04-01	400015	Antipodes Islands	Area Outside Territorial Authority	0	0	0	0	0
842606	2021-04-01	400016	Ross Dependency	Area Outside Territorial Authority	0	0	0	0	0

	SA22021_V1_00	SA22021_V1_00_NAME	SA22021_V1_00_NAME_ASCII	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	geometry	1990	1991	1992	...	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021
0	134300	Karangahape	Karangahape	0.285661	0.285661	2.602342e+03	POLYGON ((1756989.238 5919225.372, 1756983.541...	0	0	4	...	0	145	90	148	19	3	29	0	2	0
1	134500	Anzac Avenue	Anzac Avenue	0.102386	0.102386	1.499961e+03	POLYGON ((1758140.113 5920585.833, 1758170.387...	0	0	0	...	0	0	29	0	1	259	2	0	2	0
2	135100	Symonds Street North West	Symonds Street North West	0.087208	0.087208	1.610212e+03	POLYGON ((1757408.590 5919604.049, 1757394.607...	0	0	0	...	0	0	0	0	1	0	1	116	1	0
3	155400	Botany South	Botany South	0.749380	0.749380	3.657236e+03	POLYGON ((1770930.008 5910484.865, 1770924.375...	1	3	0	...	0	0	0	0	0	0	0	0	0	0
4	165600	Cape Hill	Cape Hill	0.716479	0.716479	3.571952e+03	POLYGON ((1770184.720 5882207.486, 1770188.707...	1	2	2	...	19	16	21	15	19	7	4	2	6	0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
563	170200	Pōkeno Rural	Pokeno Rural	84.945064	84.945064	7.222617e+04	POLYGON ((1784814.754 5880854.700, 1784745.051...	10	3	6	...	15	7	11	8	13	11	6	11	6	6
564	258200	Oceanic Northland Region	Oceanic Northland Region	0.000000	16289.038535	1.802429e+06	POLYGON ((1745505.932 6001808.555, 1745483.104...	0	0	0	...	0	0	0	0	0	0	0	0	0	0
565	258400	Inlet Kaipara Harbour North	Inlet Kaipara Harbour North	0.000000	330.719564	5.789059e+05	MULTIPOLYGON (((1714448.724 5982225.100, 17145...	0	0	0	...	0	0	0	0	0	0	0	0	0	0
566	258500	Oceanic Waikato Region East	Oceanic Waikato Region East	0.000000	5879.733837	8.876890e+05	POLYGON ((1804004.909 5899159.981, 1804004.056...	0	0	0	...	0	0	0	0	0	0	0	0	0	0
567	258700	Oceanic Waikato Region West	Oceanic Waikato Region West	0.000000	4184.343871	4.207090e+05	POLYGON ((1747438.207 5870452.373, 1747346.313...	0	0	0	...	0	0	0	0	0	0	0	0	0	0

	SA22021_V1_00	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	1990	1991	1992	1993	1994	1995	...	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021
count	568.000000	568.000000	568.000000	5.680000e+02	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000	...	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000	568.000000
mean	139240.322183	9.712193	76.382311	2.503863e+04	9.897887	9.198944	8.693662	10.198944	14.198944	13.464789	...	8.098592	11.165493	13.508803	16.320423	17.790493	19.274648	22.781690	26.790493	29.397887	10.297535
std	19359.762728	37.037272	821.064816	1.113122e+05	12.218187	11.223127	12.822812	17.343532	23.957580	20.031866	...	18.362338	25.395650	30.387444	37.004878	38.661058	50.248749	45.707083	49.763108	55.911654	21.721774
min	109700.000000	0.000000	0.063319	1.021312e+03	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	124175.000000	0.795242	0.824655	4.505395e+03	2.000000	2.000000	1.750000	1.000000	2.000000	2.000000	...	1.000000	1.000000	2.000000	2.000000	2.000000	3.000000	3.000000	4.000000	4.000000	0.000000
50%	138350.000000	1.124323	1.173704	5.822624e+03	7.000000	6.000000	5.000000	6.000000	7.000000	7.000000	...	3.000000	5.000000	5.000000	6.000000	6.000000	6.000000	9.000000	10.000000	13.500000	3.500000
75%	153225.000000	1.841327	2.028267	8.927031e+03	13.000000	12.000000	11.000000	12.000000	17.000000	16.000000	...	8.000000	9.000000	11.250000	12.000000	14.000000	14.250000	22.000000	29.000000	32.250000	11.000000
max	258700.000000	370.681110	16289.038535	1.802429e+06	129.000000	134.000000	145.000000	192.000000	215.000000	164.000000	...	190.000000	270.000000	288.000000	457.000000	392.000000	841.000000	674.000000	714.000000	706.000000	276.000000

	OBJECTID	STOPID	STOPCODE	STOPNAME	STOPDESC	LOCATIONTYPE	STOPLAT	STOPLON	PARENTSTATION	MODE	geometry
0	2830215	1002	1002	Stop A Lower Albert Bus Interchange	None	Stop	-36.843265	174.765754	31986	Bus	POINT (174.76575 -36.84326)
1	2830216	1003	1003	Stop B Lower Albert Bus Interchange	None	Stop	-36.843812	174.765515	31986	Bus	POINT (174.76551 -36.84381)
2	2830217	1004	1004	Stop C Lower Albert Bus Interchange	None	Stop	-36.844035	174.765507	31986	Bus	POINT (174.76551 -36.84403)
3	2830218	1005	1005	Stop D Lower Albert Bus Interchange	None	Stop	-36.843539	174.765714	31986	Bus	POINT (174.76571 -36.84354)
4	2830219	1006	1006	Stop E Lower Albert Bus Interchange	None	Stop	-36.842901	174.765329	31986	Bus	POINT (174.76533 -36.84290)
...	...	...	...	...	...	...	...	...	...	...	...
5467	2836182	8981	8981	Ngapipi Rd near Tamaki Dr	None	Stop	-36.853202	174.805741	31003	Bus	POINT (174.80574 -36.85320)
5468	2836183	8988	8988	40 Hayr Rd	None	Stop	-36.912549	174.755800	11429	Bus	POINT (174.75580 -36.91255)
5469	2836184	8990	8990	6 Carr Rd	None	Stop	-36.914679	174.755159	11430	Bus	POINT (174.75516 -36.91468)
5470	2836185	8995	8995	3 Connaught St	None	Stop	-36.926563	174.697404	11202	Bus	POINT (174.69740 -36.92656)
5471	2836186	8997	8997	21 Connaught St	None	Stop	-36.927704	174.696840	11202	Bus	POINT (174.69684 -36.92770)

	OBJECTID	STOPID	STOPCODE	STOPNAME	STOPDESC	LOCATIONTYPE	STOPLAT	STOPLON	PARENTSTATION	MODE	geometry
0	1	97	97	Papakura Train Station	None	Stop	-37.064294	174.946114	None	Train	POINT (174.94611 -37.06429)
1	2	98	98	Manurewa Train Station	None	Stop	-37.023268	174.896171	None	Train	POINT (174.89617 -37.02327)
2	3	99	99	Homai Train Station	None	Stop	-37.013463	174.874670	None	Train	POINT (174.87467 -37.01346)
3	4	100	100	Papatoetoe Train Station	None	Stop	-36.977666	174.849250	None	Train	POINT (174.84925 -36.97767)
4	5	101	101	Otahuhu Train Station	None	Stop	-36.946695	174.833210	None	Train	POINT (174.83321 -36.94669)
5	6	102	102	Penrose Train Station	None	Stop	-36.910092	174.815700	None	Train	POINT (174.81570 -36.91009)
6	7	103	103	Glen Innes Train Station	None	Stop	-36.878801	174.854120	None	Train	POINT (174.85412 -36.87880)
7	8	104	104	Morningside Train Station	None	Stop	-36.874985	174.735210	None	Train	POINT (174.73521 -36.87498)
8	9	105	105	Avondale Train Station	None	Stop	-36.897141	174.699105	None	Train	POINT (174.69910 -36.89714)
9	10	106	106	Takaanini Train Station	None	Stop	-37.041124	174.919450	None	Train	POINT (174.91945 -37.04112)
10	11	107	107	Te Mahia Train Station	None	Stop	-37.031182	174.906095	None	Train	POINT (174.90609 -37.03118)
11	12	109	109	Middlemore Train Station	None	Stop	-36.963041	174.838969	None	Train	POINT (174.83897 -36.96304)
12	13	112	112	Ellerslie Train Station	None	Stop	-36.898474	174.808111	None	Train	POINT (174.80811 -36.89847)
13	14	113	113	Greenlane Train Station	None	Stop	-36.889657	174.797419	None	Train	POINT (174.79742 -36.88966)
14	15	114	114	Remuera Train Station	None	Stop	-36.881383	174.785420	None	Train	POINT (174.78542 -36.88138)
15	16	115	115	Newmarket Train Station	None	Stop	-36.869720	174.778834	None	Train	POINT (174.77883 -36.86972)
16	17	116	116	Orakei Train Station	None	Stop	-36.862427	174.809500	None	Train	POINT (174.80950 -36.86243)
17	18	117	117	Meadowbank Train Station	None	Stop	-36.866320	174.820750	None	Train	POINT (174.82075 -36.86632)
18	19	119	119	Baldwin Ave Train Station	None	Stop	-36.877678	174.720454	None	Train	POINT (174.72045 -36.87768)
19	20	120	120	Mt Albert Train Station	None	Stop	-36.884775	174.714089	None	Train	POINT (174.71409 -36.88478)
20	21	121	121	Sunnyvale Train Station	None	Stop	-36.896868	174.632009	None	Train	POINT (174.63201 -36.89687)
21	22	122	122	Kingsland Train Station	None	Stop	-36.872508	174.744530	None	Train	POINT (174.74453 -36.87251)
22	23	123	123	Fruitvale Rd Train Station	None	Stop	-36.910670	174.667067	None	Train	POINT (174.66707 -36.91067)
23	24	124	124	Glen Eden Train Station	None	Stop	-36.910274	174.653170	None	Train	POINT (174.65317 -36.91027)
24	25	125	125	Henderson Train Station	None	Stop	-36.880965	174.630910	None	Train	POINT (174.63091 -36.88097)
25	26	126	126	Sturges Rd Train Station	None	Stop	-36.873469	174.620824	None	Train	POINT (174.62082 -36.87347)
26	27	127	127	Swanson Train Station	None	Stop	-36.866181	174.576300	None	Train	POINT (174.57630 -36.86618)
27	28	128	128	Ranui Train Station	None	Stop	-36.867937	174.603339	None	Train	POINT (174.60334 -36.86794)
28	29	129	129	New Lynn Train Station	None	Stop	-36.909385	174.684050	None	Train	POINT (174.68405 -36.90938)
29	30	130	130	Panmure Train Station	None	Stop	-36.897777	174.849670	None	Train	POINT (174.84967 -36.89778)
30	31	133	133	Britomart Train Station	None	Stop	-36.844291	174.768480	None	Train	POINT (174.76848 -36.84429)
31	32	134	134	Pukekohe Train Station	None	Stop	-37.203311	174.910160	None	Train	POINT (174.91016 -37.20331)
32	33	140	140	Parnell Train Station	None	Stop	-36.854728	174.777399	None	Train	POINT (174.77740 -36.85473)
33	34	141	141	Hamilton Frankton	None	Stop	-37.792642	175.264997	None	Train	POINT (175.26500 -37.79264)
34	35	142	142	Hamilton Rotokauri	None	Stop	-37.750795	175.229958	None	Train	POINT (175.22996 -37.75080)
35	36	143	143	Huntly	None	Stop	-37.559568	175.159498	None	Train	POINT (175.15950 -37.55957)
36	37	244	244	Sylvia Park Train Station	None	Stop	-36.914637	174.842600	None	Train	POINT (174.84260 -36.91464)
37	38	277	277	Grafton Train Station	None	Stop	-36.865612	174.769840	None	Train	POINT (174.76984 -36.86561)
38	39	605	605	Onehunga Train Station	None	Stop	-36.925394	174.786800	None	Train	POINT (174.78680 -36.92539)
39	40	606	606	Te Papapa Train Station	None	Stop	-36.920131	174.801370	None	Train	POINT (174.80137 -36.92013)
40	41	9218	9218	Manukau Train Station	None	Stop	-36.993880	174.877397	None	Train	POINT (174.87740 -36.99388)

	geometry	value
6	POINT (174.77127 -36.78666)	1
17	POINT (174.75776 -36.85657)	1
23	POINT (174.77466 -36.84825)	1
24	POINT (174.64518 -36.81060)	1
25	POINT (174.63498 -36.87765)	1
34	POINT (174.76252 -36.85307)	1
35	POINT (174.76253 -36.86491)	1
41	POINT (174.76422 -36.84391)	1
44	POINT (174.63476 -36.87823)	1
46	POINT (174.75987 -36.84590)	1
50	POINT (174.63363 -36.87871)	1
52	POINT (174.76305 -36.85285)	1
53	POINT (174.64292 -36.87171)	1
73	POINT (174.99385 -36.78114)	1
78	POINT (174.75491 -36.84145)	1
83	POINT (174.68318 -36.90920)	1
84	POINT (174.77906 -36.85296)	1
90	POINT (174.77680 -36.84937)	1
95	POINT (174.87782 -36.99102)	1
96	POINT (174.63259 -36.87919)	1
99	POINT (174.72855 -36.39045)	1
118	POINT (174.78138 -36.84959)	1
122	POINT (174.65544 -36.93813)	1
128	POINT (174.76152 -36.85897)	1
129	POINT (174.76686 -36.84898)	1
130	POINT (174.76324 -36.85088)	1

	geometry	value
0	POINT (1736573.338 5983016.712)	1
1	POINT (1748446.865 5969688.251)	1
2	POINT (1751153.511 5949091.460)	1
3	POINT (1730374.199 5939616.225)	1
4	POINT (1755337.624 5937952.716)	1
...	...	...
554	POINT (1750540.183 5930170.390)	1
555	POINT (1770436.335 5909103.133)	1
556	POINT (1770394.994 5908313.461)	1
557	POINT (1747676.110 5926691.694)	1
558	POINT (1766907.071 5897860.127)	1

	Premises_Name	Certification_Service_Type	Service_Types	Total_Beds	Premises_Website	Premises_Address_Other	Premises_Address	Premises_Address_Suburb_Road	Premises_Address_Town_City	Premises_Address_Post_Code	...	Legal_Entity_Address_Suburb_Roa	Legal_Entity_Address_Town_City	Legal_Entity_Address_Post_Code	Legal_Entity_Postal_Address	Legal_Entity_Postal_Address_Sub	Legal_Entity_Postal_Address_Tow	Legal_Entity_Postal_Address_Pos	Legal_Entity_Website	ObjectId	geometry
0	Auckland City Hospital	Public Hospital	Childrens health, Maternity, Surgical, Medical	1124	None	None	2 Park Road	Grafton	Auckland	1023	...	None	Auckland	1051	Private Bag 92024	None	Auckland	1142	http://www.adhb.govt.nz/	1	POINT (174.76948 -36.86089)
1	Auckland DHB X 3 Units - Mental Health	Public Hospital	Mental health	96	None	None	2 Park Road	Grafton	Auckland	1023	...	None	Auckland	1051	Private Bag 92024	None	Auckland	1142	http://www.adhb.govt.nz/	2	POINT (174.76948 -36.86089)
2	Buchanan Rehabilitation Centre	Public Hospital	Mental health	40	None	None	27 Sutherland Road	Point Chevalier	Auckland	1025	...	None	Auckland	1051	Private Bag 92024	None	Auckland	1142	http://www.adhb.govt.nz/	3	POINT (174.71239 -36.87277)
3	Greenlane Clinical Centre	Public Hospital	Surgical, Medical	31	None	None	214 Green Lane West	Epsom	Auckland	1051	...	None	Auckland	1051	Private Bag 92024	None	Auckland	1142	http://www.adhb.govt.nz/	4	POINT (174.77965 -36.89332)
4	Opotiki Health Care Centre	Public Hospital	Maternity, Surgical, Medical	6	None	None	32A King Street	None	Opotiki	3122	...	None	Tauranga	3112	Private Bag 12024	None	Tauranga	3143	None	5	POINT (177.28553 -38.00635)
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
78	Waitakere Hospital	Public Hospital	Geriatric, Childrens health, Surgical, Mental ...	283	None	None	55 Lincoln Road	Henderson	Auckland	0610	...	None	Auckland	622	PO Box 93503	None	Auckland	740	None	79	POINT (174.62866 -36.86952)
79	Wilson Centre	Public Hospital	Physical, Childrens health	26	None	None	1 St Leonards Road	Hauraki	Auckland	0622	...	None	Auckland	622	PO Box 93503	None	Auckland	740	None	80	POINT (174.78800 -36.79983)
80	Buller Health	Public Hospital	Maternity, Medical	8	None	None	45 Derby Street	None	Westport	7825	...	None	Greymouth	7805	PO Box 387	None	Greymouth 7840	7840	None	81	POINT (171.60491 -41.74922)
81	Grey Base Hospital	Public Hospital	Dementia care, Geriatric, Mental health, Child...	114	None	None	71 Water Walk Road	None	Greymouth	7805	...	None	Greymouth	7805	PO Box 387	None	Greymouth 7840	7840	None	82	POINT (171.19210 -42.46327)
82	Whanganui Hospital	Public Hospital	Childrens health, Maternity, Surgical, Medical...	172	None	None	100 Heads Road	Gonville	Wanganui	4501	...	None	Wanganui	4501	Private Bag 3003	None	Wanganui	4540	None	83	POINT (175.03807 -39.94539)

	name	filesize (MB)	last modified
0	bus.tif	1.26	2021-06-29 13:47:41.590
1	bus_dist.tif	2.52	2021-06-29 13:47:41.710
2	cbd.tif	1.26	2021-05-04 14:19:19.670
3	cbd_dist.tif	2.52	2021-05-04 14:19:19.720
4	cbd_dist.tifd	2.52	2021-06-29 13:47:38.110
5	consents_2018.tif	2.52	2021-06-29 15:12:05.700
6	consents_2018.tif.aux.xml	0.00	2021-06-29 15:13:06.030
7	consents_2019.tif	2.52	2021-06-29 15:12:05.710
8	consents_2020.tif	2.52	2021-06-29 15:12:05.720
9	dem.tif	5.04	2021-06-18 11:30:42.620
10	example.pkl	10.10	2021-06-29 13:47:41.720
11	hospitals.tif	1.26	2021-06-29 15:58:22.270
12	hospitals_dist.tif	2.52	2021-06-29 15:59:21.640
13	land_use_1996.tif	1.26	2021-05-04 14:45:25.500
14	land_use_2001.tif	1.26	2021-05-04 14:45:25.500
15	land_use_2008.tif	1.26	2021-05-04 14:45:25.510
16	land_use_2012.tif	1.26	2021-05-04 14:45:25.510
17	land_use_2018.tif	1.26	2021-05-04 14:45:25.520
18	mway.tif	1.26	2021-05-04 14:19:18.360
19	mway_dist.tif	2.52	2021-05-04 14:19:18.900
20	park.tif	1.26	2021-06-29 13:47:41.560
21	park_dist.tif	2.52	2021-06-29 13:47:41.470
22	pop_1996.tif	5.04	2021-05-04 14:18:13.670
23	pop_2001.tif	5.04	2021-05-04 14:18:13.680
24	pop_2006.tif	5.04	2021-05-04 14:18:13.690
25	pop_2013.tif	5.04	2021-05-04 14:18:13.700
26	pop_2018.tif	5.04	2021-05-04 14:18:13.710
27	pop_2023.tif	5.04	2021-05-04 14:18:13.720
28	pop_2028.tif	5.04	2021-05-04 14:18:13.730
29	pop_2033.tif	5.04	2021-05-04 14:18:13.740
30	pop_2038.tif	5.04	2021-05-04 14:18:13.740
31	pop_2043.tif	5.04	2021-05-04 14:18:13.760
32	pop_2048.tif	5.04	2021-05-04 14:18:13.770
33	schools.tif	1.26	2021-06-29 13:47:41.600
34	schools_dist.tif	2.52	2021-06-29 13:47:41.590
35	slope.tif	5.04	2021-06-18 11:30:42.620