GHSL#

Extract GHSL population for Sao Paulo (Brazil)#

import geopandas
import rioxarray
import contextily
import datashader as ds
import matplotlib.pyplot as plt

Functional Urban Area boundary#

The file with all the FUAs is available at:

url = ("https://jeodpp.jrc.ec.europa.eu/ftp/"\
       "jrc-opendata/GHSL/"\
       "GHS_FUA_UCDB2015_GLOBE_R2019A/V1-0/"\
       "GHS_FUA_UCDB2015_GLOBE_R2019A_54009_1K_V1_0.zip"
      )
url
'https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_FUA_UCDB2015_GLOBE_R2019A/V1-0/GHS_FUA_UCDB2015_GLOBE_R2019A_54009_1K_V1_0.zip'

From there we can build a path to read and uncompress directly on-the-fly in geopandas:

p = f"zip+{url}!GHS_FUA_UCDB2015_GLOBE_R2019A_54009_1K_V1_0.gpkg"
fuas = geopandas.read_file(p)
sao_paulo = fuas.query("eFUA_name == 'São Paulo'").to_crs("EPSG:4326")
ax = sao_paulo.plot(alpha=0.5, figsize=(9, 9))
contextily.add_basemap(ax, crs=sao_paulo.crs);
../../_images/build_ghsl_extract_7_0.png

Population layer#

The file for the tile we need is available for download at:

url = ("https://cidportal.jrc.ec.europa.eu/ftp/"\
       "jrc-opendata/GHSL/GHS_POP_MT_GLOBE_R2019A/"\
       "GHS_POP_E2015_GLOBE_R2019A_54009_250/V1-0/"\
       "tiles/"\
       "GHS_POP_E2015_GLOBE_R2019A_54009_250_V1_0_13_11.zip"
      )
url
'https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_MT_GLOBE_R2019A/GHS_POP_E2015_GLOBE_R2019A_54009_250/V1-0/tiles/GHS_POP_E2015_GLOBE_R2019A_54009_250_V1_0_13_11.zip'

From there we can build a path to read and uncompress directly on-the-fly in xarray:

%%time
p = f"zip+{url}!GHS_POP_E2015_GLOBE_R2019A_54009_250_V1_0_13_11.tif"
ghsl = rioxarray.open_rasterio(p)
ghsl
CPU times: user 40.6 ms, sys: 3.2 ms, total: 43.8 ms
Wall time: 4.56 s
<xarray.DataArray (band: 1, y: 4000, x: 4000)>
[16000000 values with dtype=float32]
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 -2e+06 -2e+06 -2.001e+06 ... -3e+06 -3e+06
  * x            (x) float64 -5.041e+06 -5.041e+06 ... -4.041e+06 -4.041e+06
    spatial_ref  int64 0
Attributes:
    _FillValue:    -200.0
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref
cvs = ds.Canvas(plot_width=600, plot_height=600)
agg = cvs.raster(ghsl.where(ghsl>0).sel(band=1))
f, ax = plt.subplots(1, figsize=(9, 7))
agg.plot.imshow(ax=ax, alpha=0.5, cmap="cividis_r")
contextily.add_basemap(
    ax, 
    crs=ghsl.rio.crs, 
    zorder=-1, 
    source=contextily.providers.CartoDB.Voyager
)
../../_images/build_ghsl_extract_13_0.png

Extraction#

To clip what we need from ghsl based on sao_paulo, we can clip it:

ghsl_sp = ghsl.rio.clip(sao_paulo.to_crs(ghsl.rio.crs).geometry.iloc[0])
ghsl_sp
<xarray.DataArray (band: 1, y: 416, x: 468)>
array([[[-200., -200., -200., ..., -200., -200., -200.],
        [-200., -200., -200., ..., -200., -200., -200.],
        [-200., -200., -200., ..., -200., -200., -200.],
        ...,
        [-200., -200., -200., ..., -200., -200., -200.],
        [-200., -200., -200., ..., -200., -200., -200.],
        [-200., -200., -200., ..., -200., -200., -200.]]], dtype=float32)
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 -2.822e+06 -2.822e+06 ... -2.926e+06 -2.926e+06
  * x            (x) float64 -4.482e+06 -4.482e+06 ... -4.365e+06 -4.365e+06
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref
    _FillValue:    -200.0

And we can write it to a GeoTIF file:

out_p = "ghsl_sao_paulo.tif"
! rm $out_p
ghsl_sp.rio.to_raster(out_p)

From rioxarray we get a COG for free!

! rio cogeo validate ghsl_sao_paulo.tif
/home/jovyan/work/data/ghsl/ghsl_sao_paulo.tif is a valid cloud optimized GeoTIFF