gNATSGO

import os
from pcxarray import pc_query, prepare_data
from pcxarray.utils import load_census_shapefile
from distributed import Client
import matplotlib.pyplot as plt

states_gdf = load_census_shapefile(level='state', verify=False)
ms_gdf = states_gdf[states_gdf['NAME'] == 'Mississippi']
ms_gdf = ms_gdf.to_crs(epsg=3814)
geom = ms_gdf.geometry.union_all()

items_gdf = pc_query(
    collections='gnatsgo-rasters',
    geometry=geom,
    crs=ms_gdf.crs,
)
print(f"Found {len(items_gdf)} items in the GNATSGO collection.")

fig, ax = plt.subplots(figsize=(5, 10))
ms_gdf.plot(ax=ax, color='lightgrey', edgecolor='black')
items_gdf.plot(ax=ax, alpha=0.2)
plt.show()

Found 10 items in the GNATSGO collection.

client = Client(processes=True)
print(f"Dask client dashboard: {client.dashboard_link}")

Dask client dashboard: http://127.0.0.1:8787/status

nccpi_da = prepare_data(
    items_gdf,
    geometry=geom,
    crs=ms_gdf.crs,
    bands=['nccpi3corn', 'nccpi3cot', 'nccpi3soy'],
    target_resolution=100, # 100 m resolution
    resampling_method='bilinear',
    chunks='auto',
    max_workers=-1,
    enable_progress_bar=True,
)
display(nccpi_da)
nccpi_da = nccpi_da.persist()

Constructing dask computation graph: 100%|██████████| 10/10 [00:00<00:00, 94.69tiles/s]

<xarray.DataArray 'reproject-1bd0dd4e259f40db9715e6cbf77d5b54' (band: 3,
                                                                y: 5393, x: 3325)> Size: 215MB
dask.array<where, shape=(3, 5393, 3325), dtype=float32, chunksize=(1, 5120, 1685), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) <U10 120B 'nccpi3corn' 'nccpi3cot' 'nccpi3soy'
  * y            (y) float64 43kB 1.578e+06 1.578e+06 ... 1.039e+06 1.039e+06
  * x            (x) float64 27kB 3.186e+05 3.186e+05 ... 6.508e+05 6.51e+05
    spatial_ref  int32 4B 3814
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0
    long_name:      National Commodity Crop Productivity Index for Corn (weig...
    _FillValue:     nan

xarray.DataArray

'reproject-1bd0dd4e259f40db9715e6cbf77d5b54'

• band: 3
• y: 5393
• x: 3325

• dask.array<chunksize=(1, 5120, 1685), meta=np.ndarray>

  Array Chunk Bytes 205.21 MiB 32.91 MiB Shape (3, 5393, 3325) (1, 5120, 1685) Dask graph 12 chunks in 131 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    <U10
    'nccpi3corn' ... 'nccpi3soy'

    array(['nccpi3corn', 'nccpi3cot', 'nccpi3soy'], dtype='<U10')

  • y
    (y)
    float64
    1.578e+06 1.578e+06 ... 1.039e+06

    units :

    metre

    resolution :

    -100.0

    crs :

    EPSG:3814

    array([1577850., 1577750., 1577650., ..., 1038850., 1038750., 1038650.],
          shape=(5393,))

  • x
    (x)
    float64
    3.186e+05 3.186e+05 ... 6.51e+05

    units :

    metre

    resolution :

    100.0

    crs :

    EPSG:3814

    array([318550., 318650., 318750., ..., 650750., 650850., 650950.],
          shape=(3325,))

  • spatial_ref
    ()
    int32
    3814

    spatial_ref :

    PROJCRS["NAD83 / Mississippi TM",BASEGEOGCRS["NAD83",DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4269]],CONVERSION["Mississippi Transverse Mercator",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",32.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-89.75,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9998335,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",1300000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["State-wide spatial data management."],AREA["United States (USA) - Mississippi."],BBOX[30.01,-91.65,35.01,-88.09]],ID["EPSG",3814]]

    crs_wkt :

    PROJCRS["NAD83 / Mississippi TM",BASEGEOGCRS["NAD83",DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4269]],CONVERSION["Mississippi Transverse Mercator",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",32.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-89.75,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9998335,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",1300000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["State-wide spatial data management."],AREA["United States (USA) - Mississippi."],BBOX[30.01,-91.65,35.01,-88.09]],ID["EPSG",3814]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    NAD83

    horizontal_datum_name :

    North American Datum 1983

    projected_crs_name :

    NAD83 / Mississippi TM

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    32.5

    longitude_of_central_meridian :

    -89.75

    false_easting :

    500000.0

    false_northing :

    1300000.0

    scale_factor_at_central_meridian :

    0.9998335

    GeoTransform :

    318500 100 0 1577900 0 -100

    array(3814, dtype=int32)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index(['nccpi3corn', 'nccpi3cot', 'nccpi3soy'], dtype='object', name='band'))

  • y
    PandasIndex

    PandasIndex(Index([1577850.0, 1577750.0, 1577650.0, 1577550.0, 1577450.0, 1577350.0,
           1577250.0, 1577150.0, 1577050.0, 1576950.0,
           ...
           1039550.0, 1039450.0, 1039350.0, 1039250.0, 1039150.0, 1039050.0,
           1038950.0, 1038850.0, 1038750.0, 1038650.0],
          dtype='float64', name='y', length=5393))

  • x
    PandasIndex

    PandasIndex(Index([318550.0, 318650.0, 318750.0, 318850.0, 318950.0, 319050.0, 319150.0,
           319250.0, 319350.0, 319450.0,
           ...
           650050.0, 650150.0, 650250.0, 650350.0, 650450.0, 650550.0, 650650.0,
           650750.0, 650850.0, 650950.0],
          dtype='float64', name='x', length=3325))

• Attributes: (5)

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

  long_name :

  National Commodity Crop Productivity Index for Corn (weighted average) for major earthy components. Values range from .01 (low productivity) to .99 (high productivity). Earthy components are those soil series or higher level taxa components that can support crop growth (Dobos et al., 2012). Major components are those soil components where the majorcompflag = 'Yes' (SSURGO component table). NULL values are presented where data are incomplete or not available.

  _FillValue :

  nan

/home/dh2306/pcxarray/.venv/lib/python3.12/site-packages/distributed/client.py:3363: UserWarning: Sending large graph of size 17.22 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

fig, ax = plt.subplots(1, 3, figsize=(15, 10), sharex=True, sharey=True)

nccpi_da.sel(band='nccpi3corn').plot.imshow(ax=ax[0], cmap='cividis', vmin=0, vmax=1, add_colorbar=False)
nccpi_da.sel(band='nccpi3cot').plot.imshow(ax=ax[1], cmap='cividis', vmin=0, vmax=1, add_colorbar=False)
nccpi_da.sel(band='nccpi3soy').plot.imshow(ax=ax[2], cmap='cividis', vmin=0, vmax=1, add_colorbar=False)

for a in ax:
    a.set_aspect('equal')
    a.set_xticks([])
    a.set_yticks([])
    a.set_xlabel('')
    a.set_ylabel('')

    # pad the axes to avoid overlap
    a.set_xlim(a.get_xlim()[0] - 5000, a.get_xlim()[1] + 5000)
    a.set_ylim(a.get_ylim()[0] - 5000, a.get_ylim()[1] + 5000)

ax[0].set_title('gNATSGO NCCPI Corn')

ax[1].set_title('gNATSGO NCCPI Cotton')

ax[2].set_title('gNATSGO NCCPI Soybean')

fig.colorbar(ax[0].images[0], ax=ax, orientation='horizontal', fraction=0.02, pad=0.02, label='National Commodity Crop Productivity Index')
plt.show()