Setup Nature Modules¶

Introduction¶

Create an example xarray.DataArray whose shape is (40 * 50).

In [1]:

import xarray as xr
import numpy as np


data = xr.DataArray(data=np.arange(2000).reshape((40, 50)), dims=("y", "x"))
data.plot()

import xarray as xr import numpy as np data = xr.DataArray(data=np.arange(2000).reshape((40, 50)), dims=("y", "x")) data.plot()

Out[1]:

<matplotlib.collections.QuadMesh at 0x157923950>

Create an ABSESpy model with 10 * 10 natural patch module.

In [ ]:

from abses import MainModel

model = MainModel()

print(f"Model created: {model}")

from abses import MainModel model = MainModel() print(f"Model created: {model}")

[15:02:11][WARNING][nature] the nature's CRS has been changed to epsg:4326.

Using ABSESpy version: v0.7.0

Out[ ]:

name     MainModel
state         init
tick             0
dtype: object

In [ ]:

from abses import load_data
from abses.viz.plotting import plot_raster

raster_path = load_data("farmland.tif")

cropland = model.nature.create_module(
    how="from_file",
    raster_file=raster_path,
    apply_raster=True,
    attr_name="test",
    name="cropland",
)

# Plot the ingested raster attribute
plot_raster(cropland, attr="test")

from abses import load_data from abses.viz.plotting import plot_raster raster_path = load_data("farmland.tif") cropland = model.nature.create_module( how="from_file", raster_file=raster_path, apply_raster=True, attr_name="test", name="cropland", ) # Plot the ingested raster attribute plot_raster(cropland, attr="test")

In [ ]:

# Create a test module for demonstration
module = model.nature.create_module(shape=(50, 40), name="test_module")

# Create a test module for demonstration module = model.nature.create_module(shape=(50, 40), name="test_module")

When apply raster of a xarray.DataArray to natural module, the raster will be automatically reprojected to the natural module's spatial extent.

In [3]:

module.apply_raster(
    data, attr_name="test", cover_crs=True, resampling_method="bilinear"
)
ax = module.get_xarray("test").plot()

module.apply_raster( data, attr_name="test", cover_crs=True, resampling_method="bilinear" ) ax = module.get_xarray("test").plot()

/Users/songshgeo/Documents/VSCode/ABSESpy/.venv/lib/python3.11/site-packages/xarray/plot/utils.py:260: RuntimeWarning: overflow encountered in scalar absolute
  vlim = max(abs(vmin - center), abs(vmax - center))

note: edge effect may exists when re-projecting the raster.

Working with real data¶

Here, we have a demo of how to load climate data and apply it to the natural module. When loading the data, the data will be automatically reprojected to the natural module's spatial extent.

In [ ]:

from abses import load_data
from abses.viz.plotting import plot_raster

raster_path = load_data("farmland.tif")

cropland2 = model.nature.create_module(
    how="from_file",
    raster_file=raster_path,
    apply_raster=True,
    name="cropland2",
    attr_name="cropland",
)

plot_raster(cropland2, attr="cropland")

from abses import load_data from abses.viz.plotting import plot_raster raster_path = load_data("farmland.tif") cropland2 = model.nature.create_module( how="from_file", raster_file=raster_path, apply_raster=True, name="cropland2", attr_name="cropland", ) plot_raster(cropland2, attr="cropland")

/Users/songshgeo/Documents/VSCode/ABSESpy/abses/nature.py:170: UserWarning: Converting PatchModule from crs EPSG:4326 to the crs of BaseNature - EPSG:4326. Please check your crs settings if this is unintended, or set `GeoSpace.warn_crs_conversion` to `False` to suppress this warning message.
  self.add_layer(module)

Out[ ]:

<matplotlib.collections.QuadMesh at 0x320f277d0>

In [5]:

data = xr.open_dataset(load_data("precipitation.nc"))
data

data = xr.open_dataset(load_data("precipitation.nc")) data

Out[5]:

<xarray.Dataset> Size: 13MB
Dimensions:  (lat: 400, lon: 700, time: 12)
Coordinates:
  * lat      (lat) float32 2kB 15.05 15.15 15.25 15.35 ... 54.75 54.85 54.95
  * lon      (lon) float32 3kB 70.05 70.15 70.25 70.35 ... 139.8 139.9 140.0
  * time     (time) datetime64[ns] 96B 2018-01-16T10:30:00 ... 2018-12-16T10:...
Data variables:
    prec     (time, lat, lon) float32 13MB ...
Attributes:
    Conventions:  CF-1.0
    history:      Created at 2016-04-02 14:55:55 UTC by the ITPCAS Data Fusio...
    source:       ITPCAS Data Fusion System (Version: 01.02.0240)
    institution:  Institute of Tibetan Plateau Research, Chinese Academy of S...
    title:        Monthly mean precipitation rate from the ITPCAS China Meteo...
    description:  Data from ITPCAS China Meteorological Forcing Dataset (Vers...
    references:   He, J., and K. Yang, 2011: China Meteorological Forcing Dat...

xarray.Dataset

• Dimensions:
  • lat: 400
  • lon: 700
  • time: 12
• Coordinates: (3)
  • lat
    (lat)
    float32
    15.05 15.15 15.25 ... 54.85 54.95

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    array([15.05    , 15.150001, 15.25    , ..., 54.75    , 54.85    , 54.95    ],
          dtype=float32)

  • lon
    (lon)
    float32
    70.05 70.15 70.25 ... 139.9 140.0

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    axis :

    X

    array([ 70.05   ,  70.15   ,  70.25   , ..., 139.75   , 139.85   , 139.95001],
          dtype=float32)

  • time
    (time)
    datetime64[ns]
    2018-01-16T10:30:00 ... 2018-12-...

    standard_name :

    time

    long_name :

    Time

    axis :

    T

    array(['2018-01-16T10:30:00.000000000', '2018-02-14T22:30:00.000000000',
           '2018-03-16T10:30:00.000000000', '2018-04-15T22:30:00.000000000',
           '2018-05-16T10:30:00.000000000', '2018-06-15T22:30:00.000000000',
           '2018-07-16T10:30:00.000000000', '2018-08-16T10:30:00.000000000',
           '2018-09-15T22:30:00.000000000', '2018-10-16T10:30:00.000000000',
           '2018-11-15T22:30:00.000000000', '2018-12-16T10:30:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • prec
    (time, lat, lon)
    float32
    ...

    long_name :

    Precipitation rate

    units :

    mm hr-1

    [3360000 values with dtype=float32]

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([15.050000190734863,  15.15000057220459,              15.25,
           15.350000381469727, 15.449999809265137, 15.550000190734863,
            15.65000057220459,              15.75, 15.850000381469727,
           15.949999809265137,
           ...
            54.04999923706055, 54.150001525878906,              54.25,
           54.349998474121094,  54.45000076293945,  54.54999923706055,
           54.650001525878906,              54.75, 54.849998474121094,
            54.95000076293945],
          dtype='float32', name='lat', length=400))

  • lon
    PandasIndex

    PandasIndex(Index([ 70.05000305175781,   70.1500015258789,              70.25,
            70.35000610351562,  70.45000457763672,  70.55000305175781,
             70.6500015258789,              70.75,  70.85000610351562,
            70.95000457763672,
           ...
            139.0500030517578, 139.14999389648438,             139.25,
           139.35000610351562, 139.45001220703125,  139.5500030517578,
           139.64999389648438,             139.75, 139.85000610351562,
           139.95001220703125],
          dtype='float32', name='lon', length=700))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-01-16 10:30:00', '2018-02-14 22:30:00',
                   '2018-03-16 10:30:00', '2018-04-15 22:30:00',
                   '2018-05-16 10:30:00', '2018-06-15 22:30:00',
                   '2018-07-16 10:30:00', '2018-08-16 10:30:00',
                   '2018-09-15 22:30:00', '2018-10-16 10:30:00',
                   '2018-11-15 22:30:00', '2018-12-16 10:30:00'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (7)

  Conventions :

  CF-1.0

  history :

  Created at 2016-04-02 14:55:55 UTC by the ITPCAS Data Fusion System

  source :

  ITPCAS Data Fusion System (Version: 01.02.0240)

  institution :

  Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China

  title :

  Monthly mean precipitation rate from the ITPCAS China Meteorological Forcing Dataset

  description :

  Data from ITPCAS China Meteorological Forcing Dataset (Version: 01.06.0014; Category: B-01: 03hr, 0.10deg)

  references :

  He, J., and K. Yang, 2011: China Meteorological Forcing Dataset. Cold and Arid Regions Science Data Center, Lanzhou, Gansu, China, digital media, doi:10.3972/westdc.002.2014.db. [Available online at http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21.]

In [6]:

a_month = data.sel(time="2013-01", method="nearest")
a_month["prec"].plot()

a_month

a_month = data.sel(time="2013-01", method="nearest") a_month["prec"].plot() a_month

Out[6]:

<xarray.Dataset> Size: 1MB
Dimensions:  (lat: 400, lon: 700)
Coordinates:
  * lat      (lat) float32 2kB 15.05 15.15 15.25 15.35 ... 54.75 54.85 54.95
  * lon      (lon) float32 3kB 70.05 70.15 70.25 70.35 ... 139.8 139.9 140.0
    time     datetime64[ns] 8B 2018-01-16T10:30:00
Data variables:
    prec     (lat, lon) float32 1MB ...
Attributes:
    Conventions:  CF-1.0
    history:      Created at 2016-04-02 14:55:55 UTC by the ITPCAS Data Fusio...
    source:       ITPCAS Data Fusion System (Version: 01.02.0240)
    institution:  Institute of Tibetan Plateau Research, Chinese Academy of S...
    title:        Monthly mean precipitation rate from the ITPCAS China Meteo...
    description:  Data from ITPCAS China Meteorological Forcing Dataset (Vers...
    references:   He, J., and K. Yang, 2011: China Meteorological Forcing Dat...

xarray.Dataset

• Dimensions:
  • lat: 400
  • lon: 700
• Coordinates: (3)
  • lat
    (lat)
    float32
    15.05 15.15 15.25 ... 54.85 54.95

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    array([15.05    , 15.150001, 15.25    , ..., 54.75    , 54.85    , 54.95    ],
          dtype=float32)

  • lon
    (lon)
    float32
    70.05 70.15 70.25 ... 139.9 140.0

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    axis :

    X

    array([ 70.05   ,  70.15   ,  70.25   , ..., 139.75   , 139.85   , 139.95001],
          dtype=float32)

  • time
    ()
    datetime64[ns]
    2018-01-16T10:30:00

    standard_name :

    time

    long_name :

    Time

    axis :

    T

    array('2018-01-16T10:30:00.000000000', dtype='datetime64[ns]')

• Data variables: (1)
  • prec
    (lat, lon)
    float32
    ...

    long_name :

    Precipitation rate

    units :

    mm hr-1

    [280000 values with dtype=float32]

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([15.050000190734863,  15.15000057220459,              15.25,
           15.350000381469727, 15.449999809265137, 15.550000190734863,
            15.65000057220459,              15.75, 15.850000381469727,
           15.949999809265137,
           ...
            54.04999923706055, 54.150001525878906,              54.25,
           54.349998474121094,  54.45000076293945,  54.54999923706055,
           54.650001525878906,              54.75, 54.849998474121094,
            54.95000076293945],
          dtype='float32', name='lat', length=400))

  • lon
    PandasIndex

    PandasIndex(Index([ 70.05000305175781,   70.1500015258789,              70.25,
            70.35000610351562,  70.45000457763672,  70.55000305175781,
             70.6500015258789,              70.75,  70.85000610351562,
            70.95000457763672,
           ...
            139.0500030517578, 139.14999389648438,             139.25,
           139.35000610351562, 139.45001220703125,  139.5500030517578,
           139.64999389648438,             139.75, 139.85000610351562,
           139.95001220703125],
          dtype='float32', name='lon', length=700))

• Attributes: (7)

  Conventions :

  CF-1.0

  history :

  Created at 2016-04-02 14:55:55 UTC by the ITPCAS Data Fusion System

  source :

  ITPCAS Data Fusion System (Version: 01.02.0240)

  institution :

  Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China

  title :

  Monthly mean precipitation rate from the ITPCAS China Meteorological Forcing Dataset

  description :

  Data from ITPCAS China Meteorological Forcing Dataset (Version: 01.06.0014; Category: B-01: 03hr, 0.10deg)

  references :

  He, J., and K. Yang, 2011: China Meteorological Forcing Dataset. Cold and Arid Regions Science Data Center, Lanzhou, Gansu, China, digital media, doi:10.3972/westdc.002.2014.db. [Available online at http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21.]

When applying the raster data to an existing module, ABSESpy automatically project match it to the current CRS.

In [7]:

cropland.apply_raster(a_month, attr_name="prec", cover_crs=True)
cropland.get_xarray("prec").plot()

cropland.apply_raster(a_month, attr_name="prec", cover_crs=True) cropland.get_xarray("prec").plot()

Out[7]:

<matplotlib.collections.QuadMesh at 0x32b8dfe10>

Dynamic raster data¶

The NetCDF dataset has three dims: 'time', 'x'(longitude), and 'y'(latitude).

One benefit of using ABSESpy is to update the dataset by real-world time. Below is an example:

In [ ]:

def update_precipitation(data, time):
    import pandas as pd
    # drop timezone to avoid mismatch
    idx = pd.DatetimeIndex(data.indexes["time"]).tz_localize(None)
    ds = data.assign_coords(time=idx)
    # use a naive datetime for selection
    return ds.sel(time=time.dt.naive(), method="nearest")


cropland.add_dynamic_variable(
    name="prec",
    data=data,
    function=update_precipitation,
    cover_crs=True,
)

def update_precipitation(data, time): import pandas as pd # drop timezone to avoid mismatch idx = pd.DatetimeIndex(data.indexes["time"]).tz_localize(None) ds = data.assign_coords(time=idx) # use a naive datetime for selection return ds.sel(time=time.dt.naive(), method="nearest") cropland.add_dynamic_variable( name="prec", data=data, function=update_precipitation, cover_crs=True, )

In [9]:

cropland.dynamic_var("prec", dtype="xarray").plot()

cropland.dynamic_var("prec", dtype="xarray").plot()

Out[9]:

<matplotlib.collections.QuadMesh at 0x3233414d0>

Create Module from Vector Dataset¶

In this example, we will create a PatchModule and geo-agents with a shapefile dataset.

In [10]:

import geopandas as gpd

data_path = load_data("YR_cities.zip")
gdf = gpd.read_file(data_path)
gdf.head()

import geopandas as gpd data_path = load_data("YR_cities.zip") gdf = gpd.read_file(data_path) gdf.head()

Out[10]:

	City_ID	Ratio	area	Shaoefile_	Perfecture	Province_n	geometry
0	100	0.546595	3958.281245	128.0	C100	Henan	POLYGON ((113.38580 35.47040, 113.39200 35.464...
1	102	0.385113	15258.027194	170.0	C102	Henan	POLYGON ((112.03220 35.04700, 112.03260 35.045...
2	106	0.567869	4273.313471	197.0	C106	Henan	POLYGON ((115.19540 36.20580, 115.19700 36.203...
3	107	0.389969	9951.746614	215.0	C107	Henan	POLYGON ((111.96140 35.07840, 111.96400 35.076...
4	109	0.601268	8279.739196	290.0	C109	Henan	POLYGON ((113.67120 35.83740, 113.67660 35.835...

In [ ]:

# Create a new model, and import the above GeoDataFrame as a patch module.
model = MainModel()
yr_basin = model.nature.create_module(
    how="from_vector",
    vector_file=gdf,
    resolution=0.1,  # 0.1 degree.
    major_layer=True,
    name="yellow_river",
    attr_name="area",
    apply_raster=True,
)
from abses.viz.plotting import plot_raster
plot_raster(yr_basin, attr="area")

# Create a new model, and import the above GeoDataFrame as a patch module. model = MainModel() yr_basin = model.nature.create_module( how="from_vector", vector_file=gdf, resolution=0.1, # 0.1 degree. major_layer=True, name="yellow_river", attr_name="area", apply_raster=True, ) from abses.viz.plotting import plot_raster plot_raster(yr_basin, attr="area")

[15:02:13][WARNING][nature] the nature's CRS has been changed to GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]].
/Users/songshgeo/Documents/VSCode/ABSESpy/abses/nature.py:170: UserWarning: Converting PatchModule from crs EPSG:4326 to the crs of BaseNature - EPSG:4326. Please check your crs settings if this is unintended, or set `GeoSpace.warn_crs_conversion` to `False` to suppress this warning message.
  self.add_layer(module)
/Users/songshgeo/Documents/VSCode/ABSESpy/abses/viz/viz_nature.py:56: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax.legend(**legend_kwargs)

Out[ ]:

<Axes: title={'center': 'spatial_ref = 0'}, xlabel='x', ylabel='y'>

In [ ]:

def update_precipitation(data, time):
    import pandas as pd
    idx = pd.DatetimeIndex(data.indexes["time"]).tz_localize(None)
    ds = data.assign_coords(time=idx)
    return ds.sel(time=time.dt.naive(), method="nearest")


yr_basin.add_dynamic_variable(
    name="prec",
    data=data,
    function=update_precipitation,
    cover_crs=True,
)

def update_precipitation(data, time): import pandas as pd idx = pd.DatetimeIndex(data.indexes["time"]).tz_localize(None) ds = data.assign_coords(time=idx) return ds.sel(time=time.dt.naive(), method="nearest") yr_basin.add_dynamic_variable( name="prec", data=data, function=update_precipitation, cover_crs=True, )

In [13]:

from abses import Actor


class County(Actor):
    """City in the Yellow River Basin."""

    marker = "^"
    color = "#FFFF80"


model.agents.new_from_gdf(
    gdf,
    unique_id="City_ID",
    attrs=["area", "Province_n"],
)
yr_basin.random.new(County, size=50)

from abses import Actor class County(Actor): """City in the Yellow River Basin.""" marker = "^" color = "#FFFF80" model.agents.new_from_gdf( gdf, unique_id="City_ID", attrs=["area", "Province_n"], ) yr_basin.random.new(County, size=50)

Out[13]:

<ActorsList: (50)County>

In [ ]:

from abses.viz.plotting import plot_raster
plot_raster(yr_basin, attr="prec")

from abses.viz.plotting import plot_raster plot_raster(yr_basin, attr="prec")

Out[ ]:

<Axes: title={'center': 'spatial_ref = 0'}, xlabel='x', ylabel='y'>