Satellite-Guided Delineation of Crop Production Zones from Official Crop Statistics for Spatial Agricultural Decision Support

1. Introduction

Spatial information on crop yields is central for assessing climate risks, evaluating management strategies, and designing policies for water and nutrient management. Yet in most regions, the only consistent, long-term production statistics are reported at administrative units such as counties or districts. In the United States, the USDA National Agricultural Statistics Service (NASS) provides annual county-level planted area, harvested area, and yield for major crops, which underpin a huge body of agronomic and economic work [1]. However, these data say nothing about how production is distributed within counties, or how irrigated “cores” differ from marginal dryland areas scales that matter for both process-based models and on-farm decision support.

Remote sensing has enabled substantial progress in estimating crop yields below administrative reporting units, particularly through vegetation indices, meteorological variables, crop models, and statistical or machine-learning approaches. These studies have demonstrated that satellite-derived indicators can capture important spatial and temporal variation in crop growth and productivity, especially when calibrated against plot measurements, combine-harvester yield monitors, or well-parameterized crop models [2,3,4,5,6,7]. For example, the Scalable Satellite-based Crop Yield Mapper (SCYM) framework combines crop-model simulations and Landsat-scale vegetation indices to estimate maize and soybean yields across the U.S. Corn Belt [3], while field-scale studies have fused satellite imagery with combine-yield observations to generate 10–30 m yield maps for wheat and other cereals [8,9,10]. These approaches are valuable for yield prediction and precision-agriculture applications, but they address a different problem from the one considered here. They generally estimate yield independently from remote-sensing or model-derived predictors and then compare predictions with official statistics or field observations. In contrast, many regional modeling and decision-support applications require spatial yield layers that remain explicitly consistent with official harvested area and county-level yield statistics, particularly where field-level yield observations are sparse, proprietary, or unavailable.

Beyond purely statistical remote-sensing approaches, a large body of work relies on process-based crop growth models such as WOFOST, APSIM, and DSSAT, sometimes coupled with data assimilation of satellite greenness indices or soil moisture to forecast yields from local to continental scales [11,12,13,14]. These models explicitly represent soil–plant–atmosphere processes and management, but their operational use is often limited by the need for detailed input data and careful calibration in each new region. In parallel to process-based models, many studies have pursued purely data-driven routes in which yields are inferred directly from satellite time series and ancillary weather data. Some approaches rely on vegetation indices at a few key phenological stages or on phenology metrics derived from VI curves to forecast county-scale yields weeks to months before harvest [15,16]. Others exploit full-season Landsat or Sentinel-2 trajectories, sometimes combined with gridded climate variables, to calibrate empirical or machine-learning models at field scale [17]. Reviews of these machine-learning approaches emphasize the flexibility of algorithms such as random forests, gradient boosting, and Gaussian processes to handle multi-sensor inputs and nonlinear responses, while also noting their dependence on large, high-quality yield datasets for training and validation [18,19]. Despite their successes, these model- and ML-based approaches generally require either dense local yield observations or rich management datasets, and they typically use official statistics (e.g., NASS) only as a loose consistency check rather than as a hard constraint on total production.

A parallel stream of work uses remote sensing to construct gridded yield or production surfaces at regional to global scales by downscaling administrative statistics. Classic examples include Spatial Production Allocation Model (SPAM) [20] and M3-Crops [21], which distribute national or subnational crop statistics to grid cells based on cropland maps and suitability layers. More recently, Qin, Wu, Zeng, Zhang and Tian [5] developed the GGCP10 dataset, a global gridded crop-production product at 10 km resolution for 2010–2020 that integrates multiple satellite, climate, and soil datasets with machine-learning models, and then calibrates grid-level production back to FAO national statistics. Complementing these statistical products, the Global Gridded Crop Model Intercomparison (GGCMI) and related AgMIP studies use ensembles of process-based models (e.g., DSSAT, LPJmL, PRYSBL2) to simulate yields on 0.5° grids under historical and future climates, providing multi-model projections of yield responses to warming and CO₂ [22]. These products demonstrate the feasibility of enforcing consistency with official statistics while exploiting remote-sensing signals, but they operate at relatively coarse spatial resolutions and often aggregate management environments (e.g., irrigated vs. rainfed) within a grid cell.

Despite these advances, a methodological gap remains between two dominant types of yield information. On one side, official agricultural statistics such as NASS provide trusted, consistent, long-term estimates of crop area and yield, but only at coarse administrative scales. On the other side, field- or pixel-scale yield-mapping approaches can resolve fine spatial variability, but they often require local yield observations, proprietary yield-monitor data, detailed management information, or crop-model calibration. Existing gridded crop-production datasets and downscaling products demonstrate the value of spatializing agricultural statistics, but they typically operate at coarse spatial resolutions or rely on global suitability, climate, or machine-learning layers that are not designed to delineate local production zones within individual counties. What is still missing is a transparent, reproducible, and constraint-preserving workflow that uses public satellite data to redistribute official county-level statistics into spatially explicit production units while preserving both harvested area and county mean yield by construction. Furthermore, spatial disaggregation has become an important geospatial modeling problem when detailed spatial processes must be inferred from aggregated administrative statistics, particularly when large spatial domains and fine target resolutions make computation challenging [23].

At the same time, the U.S. now has a mature set of wall-to-wall land-use and remote-sensing products that, in principle, could support systematic disaggregation of NASS yields. The Cropland Data Layer (CDL) provides annual crop maps at 30 m resolution for the conterminous U.S. and has been widely used to delineate crop-specific masks for yield-gap and production analyses. Multi-sensor NDVI products from MODIS, Landsat, and Sentinel-2 offer dense time series of canopy greenness, while satellite-based evapotranspiration (ET) products such as OpenET SSEBop CONUS GridMET v2.0 quantify spatial patterns in water use. Together with the NASS Quick Stats API, which exposes county statistics programmatically via tools such as the rnassqs R package, these datasets create an opportunity: instead of fitting yield models to ground-truth data, we can treat the county mean as fixed and use remote sensing only to decide how that mean should be distributed across space.

This paper develops and demonstrates nayd R Package, a NASS-anchored yield-disaggregation workflow that converts county-level harvested area and yield statistics into spatially explicit production units and zonal clusters using only publicly available datasets. The goal is not to estimate absolute yield independently of NASS, but to spatially redistribute official crop statistics in a way that preserves both harvested area and county mean yield by construction. The workflow: (i) refines CDL crop masks using an NDVI-based eligibility filter so that the retained crop area is consistent with NASS harvested area; (ii) combines seasonal NDVI and satellite-derived evapotranspiration into normalized spatial weights representing relative within-county production potential; and (iii) aggregates weighted pixels into realistic harvest units and graph-based production clusters, then allocates county-level yield to these units under explicit area and yield constraints. All steps are implemented in open-source R code and are designed to be adaptable across crops and regions where crop masks, aggregated yield statistics, vegetation-index products, and satellite-derived water-use indicators are available.

We illustrate the workflow for cotton and winter wheat in two Texas county-year case studies and evaluate its external plausibility using independent Texas A&M AgriLife variety-trial data across multiple years and locations. The workflow is not fitted to the trial data; therefore, the comparison is used as an independent check of whether NASS-anchored production clusters reproduce broad spatial yield patterns rather than exact field-level performance.

The specific objectives were to: (1) develop and document an open-source, constraint-preserving workflow for converting county-level crop statistics into spatially explicit production units and zonal clusters; (2) evaluate how CDL crop masks, NDVI-based eligibility filtering, and NDVI–ET weighting determine effective harvested area and within-county yield gradients; and (3) assess whether the resulting cluster yields agree with independent cotton and winter wheat variety-trial data without using those trials for calibration. The main contribution is a transparent geospatial allocation framework that bridges official agricultural statistics and satellite-derived spatial variability while preserving aggregate statistical consistency.

2. Materials and Methods

2.1. Study Region

This study illustrates cotton and winter wheat NASS–yield disaggregation in two counties in the Texas High Plains (THP), a semi-arid agricultural region in northwestern Texas covering approximately 8.9 million ha. The THP is characterized by a relatively flat terrain and a continental climate with hot summers, cold winters, and erratic precipitation of roughly 450–550 mm annually. Major crops include cotton, maize, sorghum (Sorghum bicolor L. Moench), and winter wheat, with cotton relying heavily on irrigation and wheat more often grown under rainfed conditions in the northern and central counties [24]. Ongoing Ogalala aquifer depletion has driven a gradual shift toward more water-efficient and rainfed systems, increasing the region’s vulnerability to climate variability and water constraints.

Within this broader region, the workflow is illustrated for cotton in Gaines County in 2021 and for winter wheat in Deaf Smith County in 2016. To evaluate the plausibility of the disaggregated yields beyond these examples, we additionally compared cluster mean yields with variety-trial data for cotton (2018–2024) and winter wheat (2008–2023) from Texas A&M AgriLife programs at a subset of locations across Texas where trials were available. These sites span multiple counties and environments, including both irrigated and rainfed systems, and thus provide an independent check on whether the NASS-anchored clusters reproduce broad spatial patterns in yield.

2.2. Data Sources

2.2.1. NASS County Statistics

County-level planted area, harvested area, and cotton lint yield and winter wheat grain yield were obtained from the USDA National Agricultural Statistics Service (NASS) for Texas counties in the needed years using the rnassqs R package, which provides an interface to the NASS ‘Quick Stats’ web API (https://quickstats.nass.usda.gov/api/) with built-in functions that facilitate building queries based on available parameters. Planted and harvested areas (ha) were used as constraints on the spatial masks (Section 2.3.1), and county average yield (kg ha⁻¹) provided the target for disaggregation to field and cluster scales (Section 2.3.4 and Section 2.3.5.

2.2.2. Cropland Data Layer (CDL)

Annual USDA Cropland Data Layer (CDL) maps were used to identify crop-specific pixels for cotton and winter wheat from 2008 onwards (https://www.nass.usda.gov/Research_and_Science/Cropland/Release/index.php). CDL rasters were reprojected to an equal-area projection (NAD83 / Conus Albers; EPSG:5070) and resampled to a common grid. Pixels mapped as the target crop in the CDL were combined into a binary crop mask with an optional user-defined lower and upper buffer relative to NASS planted area. These buffers are specified as multiplicative factors on area (e.g., 0.6 and 2.3, corresponding to 60% and 230% of the NASS planted area) to accommodate CDL underestimation and overestimation while excluding obvious misclassified patches. When CDL underestimation exceeds a user-defined threshold (e.g., >40%), the workflow expands the mask by taking the union of crop pixels across the three years surrounding the target year; when CDL overestimation exceeds a user-defined threshold (e.g., >70%), it tightens the mask by taking the intersection of crop pixels across the three surrounding years. In our Texas cases, CDL underestimation relative to NASS planted area was uncommon, while overestimation occurred more frequently; nonetheless, setting the upper buffer too low can eliminate genuine crop pixels and should be avoided.

2.2.3. Remote Sensing: NDVI and ET

Multi-sensor vegetation indices and satellite-based evapotranspiration (ET) were used to construct spatial weights for yield disaggregation. For vegetation, we used NDVI derived from MODIS 16-day collections "MODIS/061/MYD13Q1" (250 m, 2023–present) or "MODIS/061/MOD13Q1" (250 m, 2000–2022), or Landsat Tier 1 Collection 2 surface reflectance from Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI (30 m, 1983–present), or Sentinel-2 Level-2A surface reflectance from Sentinel-2A/B MSI (10 m, 2016–present), depending on sensor availability and user preference for each county–year. All products were accessed as atmospherically corrected surface reflectance. Quality control was performed using the accompanying QA bands; (i) MODIS MO/YD13Q1: we excluded pixels flagged as clouds, cloud shadows, snow/ice, or other low-quality observations in the Summary QA and pixel reliability flags; (ii) Landsat: we used the pixel_qa band from Collection 2 Tier 1 SR to remove water, clouds, cloud shadows, and snow; (iii) Sentinel-2: we masked clouds and cloud shadows using the Scene Classification (SCL) band and/or quality bits, e.g., classes “cloud”, “cloud shadow”, “cirrus”, and “snow/ice”. Only cloud-free, snow-free observations over land were retained for NDVI compositing.

NDVI was calculated from red and near-infrared reflectance as:

$$NDVI={\displaystyle \frac{\mathsf{\rho }\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{ }\mathsf{\rho }\mathrm{R}\mathrm{E}\mathrm{D}}{\mathsf{\rho }\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{ }\mathsf{\rho }\mathrm{R}\mathrm{E}\mathrm{D}\mathrm{ }}}$$

(1)

Within the workflow, the Spatial_Weights module constructs the NDVI–ET weight surface and allows the user to select which NDVI sensor (MODIS, Landsat, or Sentinel-2) to use in a given county-year. This selection is subject to the sensor-specific availability and quality screening (e.g., MODIS from 2000 onward, Sentinel-2 from 2016 onward, and only cloud- and snow-free observations retained). The NDVI composites were generated over crop-specific growing-season windows and reprojected to the common 30 m grid in EPSG:5070 using bilinear interpolation when downscaling from coarser products, e.g., MODIS 250 m, or upscaling the original Sentinel-2 scenes. Season start and end dates are user-defined. In this case study, we set the NDVI compositing windows from 1 June to 30 September for cotton and from 1 March to 15 May for winter wheat. These windows ensure that the target crop has emerged and reached canopy development sufficient to cover the soil. In THP, cotton is planted by mid-May and winter wheat resumes spring growth by March.

For ET, we used monthly actual ET from the OpenET (SSEBop CONUS GRIDMET v2.0) Landsat-based 30 m spatial resolution product and aggregated to the seasonal totals over the same crop-specific windows for each crop and year. All satellite products were resampled and reprojected to a common 30 m grid in EPSG:5070, using nearest-neighbor resampling for categorical layers and bilinear resampling for continuous layers (NDVI, ET). Within each season, mean NDVI and summed ET were computed for each pixel to represent crop vigor and water use. These two signals were then combined into a single weight layer using user-defined mixture coefficients, e.g., 70% weight on NDVI and 30% weight on ET, as described in the following sections.

2.3. Software Workflow

All components of the NASS-anchored yield disaggregation workflow are implemented in open-source R code and provided in a public GitHub repository as an R library nayd (Table 1). In what follows, we refer to this implementation simply as “the workflow”. Table 2 describes the required libraries and their sources for the nayd library.

County-level NASS statistics are combined with publicly available remote-sensing products to generate spatially explicit production units and zonal clusters (Figure 1).

First, CDL crop masks are refined using an NDVI-based eligibility filter so that the retained crop area matches NASS harvested area within a user-defined tolerance. Second, seasonal NDVI and satellite-derived evapotranspiration (ET) are normalized and combined into spatial weights representing the relative contribution of each pixel to the county mean yield. Third, these weights are aggregated into realistic harvest units and grouped into zonal production clusters using a graph-based connectivity rule with a user-defined distance threshold (10 km in this study). Finally, NASS county mean yield is disaggregated to units and clusters using weighted allocation with a tunable contrast parameter (β), producing spatially explicit yield surfaces suitable for regional modeling and validation.

2.3.1. Crop Masks and NDVI Eligibility

In the workflow, the USDA CDL is used as the initial crop mask for the target crop from 2008 onward. Building on this CDL mask, we implemented an NDVI eligibility function to further restrict the mask to pixels with an actual crop signal (Section 2.2.2). The purpose of this step is to tune the effective crop area so that it is consistent with the NASS harvested area.For each county-year and crop, we applied the NDVI eligibility filter as follows. We first computed a seasonal NDVI composite over pixels classified as the target crop in the CDL. We then search over a grid of NDVI quantiles $p$(from high to low, e.g. 0.90, 0.88, …, 0.10). For each candidate quantile $p$, the function: (i) computes an NDVI threshold $T_p$as the $p$-th quantile of seasonal NDVI over CDL crop pixels; (ii) retains pixels with NDVI $\ge T_p$; (iii) computes the total area of retained pixels using per-cell areas (via terra::cellSize in EPSG:5070); and (iv) compares this area to a target $A_{\mathrm{target}}=A_{\mathrm{harv}}\times (1+f_{\mathrm{buf}})$, where $A_{\mathrm{harv}}$is the NASS harvested area and $f_{\mathrm{buf}}$is a user-defined fractional buffer (e.g. 0.10-0.20, corresponding to 10-20%). The procedure stops at the first quantile for which the retained area exceeds $A_{\mathrm{target}}$.

The resulting NDVI eligibility mask removes low-signal or misclassified CDL pixels and produces an effective crop area that closely approximates the NASS harvested area (within the chosen buffer). This mask defines the spatial support for all subsequent weighting and disaggregation steps in the workflow.

2.3.2. Construction of NDVI-ET Weight Layer

For each county-year and crop, a spatial weight surface was constructed combining NDVI and ET information over the eligible pixels. First, NDVI and ET layers were separately min-max scaled to a common weight range. For each pixel i in county c and year y,

$${\mathit{w}}_{i,c,y\mathrm{ }}^{NDVI}={\displaystyle \frac{NDV{I}_{i,c,y}-\mathrm{m}\mathrm{i}\mathrm{n}\left(NDV{I}_{c,y}\right)}{\max \left(NDV{I}_{c,y}\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(NDV{I}_{c,y}\right)}}\times \left(1-\mathrm{ }\in \right)+\in$$

(2)

$${\mathit{w}}_{i,c,y}^{ET}={\displaystyle \frac{E{T}_{i,c,y}-\min \left(E{T}_{c,y}\right)}{\max \left(E{T}_{c,y}\right)-\min \left(E{T}_{c,y}\right)}}\times \left(1-\mathrm{ }\in \right)+\in$$

(3)

where $w_{i,c,y}^{NDVI}and{w}_{i,c,y}^{ET}$are the normalized NDVI and ET weights, respectively, and $\in$ is a small positive constant, e.g. 0.01, that prevents 0 weights. The normalized NDVI and ET weights layers were then combined according to a user-defined mixture parameters, ${\alpha }_{NDVI}$ and ${\alpha }_{ET}$, therefore the raw weight equals

$${\mathit{w}}_{raw}={\alpha }_{NDVI}\times {\mathit{w}}_{NDVI}+{\alpha }_{ET}\times {\mathit{w}}_{ET}$$

(4)

where $w_{raw}$ is the raw weight of NDVI and ET mixed layer, and ${\alpha }_{NDVI}and{\alpha }_{ET}$ are user-defined parameters for NDVI and ET, respectively (e.g., 70% for NDVI and 30% for ET). The sum of ${\alpha }_{NDVI}and{\alpha }_{ET}$ must equal to 1, but if the ET is unavailable, i.e. < 1999, ${\mathit{w}}_{raw}$ is set to ${\mathit{w}}_{NDVI}$. Finally, the ${\mathit{w}}_{raw}$ was renormalized ${\mathit{w}}_{norm}$so that the sum of weights over eligible pixels would be equal to 1. The resulting final normalized weight raster represents a continuous probability surface over the NDVI eligible crop area, with higher values in pixels more likely to contribute to the county’s harvested yield.

2.3.3. Segmentation Into Production Units

Within each county-year and crop, the normalized weight ${\mathit{w}}_{norm}$was aggregated into spatial production units using a simple segmentation and tiling procedure. First, contiguous patches of NDVI-eligible pixels were delineated using the patches function in the terra package with 4-neighbor connectivity. Tiny specks (< 50 cells) were removed with a sieve operation to ensure that remaining patches represented meaningful contiguous crop areas. For each patch, we computed its total area (ha) using cellSize in EPSG:5070 and the sum of normalized weights from ${\mathit{w}}_{\mathrm{norm}}$.

Patches were then divided into two classes based on a user-defined “giant field” threshold ${\mathit{A}}_{\mathrm{giant}}$(ha). Patches with area $\le {\mathit{A}}_{\mathrm{giant}}$were treated as standalone production units, whereas patches with area $>{\mathit{A}}_{\mathrm{giant}}$were further subdivided into pseudo-fields using a regular fishnet. In this case study we set ${\mathit{A}}_{\mathrm{giant}}\approx 800$ha for both cotton and winter wheat, reflecting the presence of very large, consolidated fields in parts of the THP, but this parameter is fully tunable in the workflow.

For “giant” patches, we generated a square fishnet with a target tile area ${\mathit{A}}_{\mathrm{tile}}$(default ≈ 200 ha), clipped the grid to the patch boundary, and for each resulting tile computed area (ha), total weight, and weight density. Each tile was then treated as a separate pseudo-field nested within the original patch. Finally, standalone patches and tiles were merged into a single set of production units and filtered to retain units within a user-specified size range, e.g. 5-1000 ha in our case study. The remaining units were sorted in descending order of total weight and assigned a cumulative area, which forms the basis for selecting harvest units that match the NASS harvested area.

Small and medium fields and all tiles were then combined into a single set of production units. For each unit $u$, we computed its area ${\mathit{A}}_u$(ha), total weight ${\mathit{W}}_u$(the sum of ${\mathit{w}}_{\mathrm{norm}}$over all pixels in the unit), and weight density ${\mathit{W}}_u/{\mathit{A}}_u$. Units were filtered to retain only those within a plausible size range specified by the user. In this case study we kept units between 5 and 1000 ha, which removes very small slivers and extremely large polygons. These lower and upper bounds are exposed as arguments in the workflow so that users can adapt them to other crops and regions. The remaining units were sorted in descending order of total weight ${\mathit{W}}_u$and assigned a cumulative area:

$$cum\_are{a}_k={\sum }_{u=1}^k{\mathit{A}}_{\mathrm{u}}$$

(5)

which forms the basis for selecting the subset of “harvest units” whose total area matches the NASS harvested area.

2.3.4. Production Units Subset and NASS-Anchored Disaggregation

Within each county-year, a subset of production units (u) was selected such that their cumulative area approximated the A_harv, by ranking ${\mathit{W}}_u$ and calculating the corresponding $cu{m}_{are{a}_k}$and identifying the index at which the $cu{m}_{are{a}_k}$first exceeded the A_harv. The selected units were rasterized to a mask and used to restrict the ${\mathit{w}}_{norm}$ layer to A_harv. The restricted weights were then renormalized so that their sum would be equal to 1, and the pre-unit weight sums over the harvest subset were recomputed:

$$A_u^{*}=A_u\times {\displaystyle \frac{A_{harv}}{{\mathsf{\Sigma }}_{\mathrm{u}}{A}_u}}$$

(6)

where $A_u^{*}$ is the scaled analytical area of unit $u$. This scaling preserves the spatial geometry of the selected production units while ensuring that the total analytical support is consistent with the official harvested area. Therefore, the disaggregation is constrained by both the NASS harvested area and the NASS county mean yield.

To control the spread of yields around the NASS mean while preserving the NASS anchor, a shape parameter was added $\left(\beta >=0\right)$. For each unit u in the harvest subset, a transformed weight was computed

$$r_u=W_u^{\beta }$$

(7)

where $r_u$is the beta-transformed relative productivity score for unit $u$. When $\beta =1$, the yield pattern follows the original NDVI–ET-derived unit weights. Values of $\beta >1$ amplify spatial contrasts by increasing the relative influence of high-weight units, whereas values of $\beta <1$compress the distribution and produce a flatter yield pattern across units. In the experiments below, several $\beta$values were evaluated to quantify how this parameter affects agreement with independent variety-trial data.

Importantly, the transformed scores were not treated as direct yield values or simple production shares. Instead, they were interpreted as relative productivity multipliers and normalized by their area-weighted mean:

$$\overline{r_A}={\displaystyle \frac{{\Sigma }_u{r}_u{A}_u^{*}}{{\Sigma }_u{A}_u^{*}}}$$

(8)

and unit-level raw yield was then calculated as:

$$\mathrm{ }{\mathrm{Y}}_{\mathrm{u}}^{\mathrm{r}\mathrm{a}\mathrm{w}}=Y_{NASS}\mathrm{ }\left({\displaystyle \frac{r_u}{r_A}}\right)$$

(9)

where $Y_{NASS}$ is the NASS-reported county mean yield. This formulation guarantees that, before applying yield safeguards, the area-weighted mean of unit-level yields equals the NASS county mean:

$${\displaystyle \frac{{\sum }_u{Y}_u^{raw}{A}_u^{*}}{{\sum }_u{A}_u^{*}}}=Y_{NASS}$$

(10)

Thus, the NASS mean yield is preserved by construction while allowing spatial variation among production units according to the NDVI–ET-derived relative productivity scores. To avoid unrealistically low or extreme yield values, crop-specific lower and upper yield bounds were applied to the raw unit-level yields:

$$Y_u^{clamped}=\mathrm{m}\mathrm{i}\mathrm{n}[\max \left(Y_u^{raw},\mathrm{ }{Y}_{floor}\right),\mathrm{ }{Y}_{ceiling}]$$

(11)

where $Y_{floor}$ and $Y_{ceiling}$ represent crop-specific and, when used, relative yield limits. Because this clipping step can slightly shift the area-weighted mean away from $Y_{NASS}$, the clamped yields were re-centered using an adjustment factor:

$$f_{adj}={\displaystyle \frac{Y_{NASS}}{\frac{{\sum }_{u\mathrm{ }}{Y}_u^{clamped\mathrm{ }}{A}_u^{\mathrm{*}}}{{\sum }_u{A}_u^{\mathrm{*}}}}}$$

(12)

and final unit-level yields were calculated as:

$$Y_u^{final}=Y_u^{clamped}\mathrm{ }\times f_{adj}\mathrm{ }$$

(13)

A final clipping step was applied when necessary to ensure that adjusted yields remained within the specified crop-specific bounds. Unit-level production was then calculated as:

$$P_u=Y_u^{final}{A}_u^{*}$$

(14)

The resulting production-unit yields and production totals therefore preserve the NASS harvested-area and county-yield constraints as closely as possible while expressing within-county spatial variability derived from NDVI and ET. Any residual deviation from the NASS county mean after clipping and re-centering was reported in the diagnostic summaries.

2.3.5. Construction of Production Clusters

To generate coarser production zones suitable for coupling with crop models or gridded soil datasets, the production units (u) were grouped into production clusters using a spatial connectivity rule. The centroids of harvest units were treated as nodes in a graph. Two units $u_i$ and $u_j$were connected if the Euclidean distance between their centroids, d($u_i,u_j),$was <= a user-defined clustering radius d_clust. In this study, we set $d_{\mathrm{clust}}=10$km. Clusters were then defined as the connected components of this graph, so each production cluster consists of harvest units linked by chains of neighbors separated by at most $d_{\mathrm{clust}}$. Individual cluster areas vary depending on field density. The 10 km was chosen to suit large scale grid data resolution. Nonetheless, end-users have the option to either decrease or increase the clustering distance to suit their objectives.

For each cluster $c$, the areas and weights of its member units were summed as follows:

$${\mathit{A}}_c= \sum _{u\in c}{\mathit{A}}_u$$

(15)

$${\mathit{W}}_c=\sum _{u\in c}{\mathit{W}}_u$$

(16)

Cluster-level yields were obtained by aggregating unit yield Y_u within each cluster, i.e., area-weighted mean, ensuring that the area-weighted mean of cluster yields equals the NASS county-level yield. The resulting cluster yields represent production zones at approximately 10 km spatial scale.

2.4. Evaluation Against Variety Trials

We evaluated the plausibility of disaggregated yields using independent variety trial data for cotton and winter wheat obtained from the Texas A&M AgriLife Variety Testing Programs. For cotton, we used irrigated and dryland Realistic Agronomic Cotton Evaluation (RACE) trials in Texas from 2018 to 2024 (https://varietytesting.tamu.edu/cotton/), extracting the test mean lint yield for each site-year. For winter wheat, we used the Texas Small Grains variety trial reports (https://varietytesting.tamu.edu/smallgrains/) from 2008 to 2023 and extracted the site mean grain yield. Each site was characterized by year, county, environment (irrigated or dryland), and either exact coordinates or the nearest town/experiment station, from which coordinates were approximated.

After extracting the variety trial data and running the workflow for the corresponding county-year, each trial was matched to the disaggregated cluster yields within search radii of 10, 20, and 30 km, selecting the cluster whose mean yield was closest to the observed trial mean. These radii are appropriate because we do not validate individual fields but broader production clusters: each cluster aggregates many fields over several kilometers. In practice, this means that for each trial point we ask whether the clustering produces any nearby production zone, within a given radius, whose mean yield is consistent with the trial yield. The search was restricted to clusters from the same county and year; distances between the trial coordinates and cluster centroids were computed in EPSG:5070, candidate clusters within each search radius were identified, and among these candidates we selected the cluster minimizing the absolute difference between cluster yield and trial yield. These radii reflect both the typical spatial scale of the production clusters and location uncertainty in the trial coordinates.

In addition, several β values (eq. 7) were tested (β = 0.8, 1.0, 1.2, …, 1.8) to determine which combination of contrast parameter and search radius provided the best agreement with the variety trial data. Thus, for each crop we evaluated 18 configurations (6 β values × 3 radii). Comparisons between variety trial yields and cluster-level disaggregated yields were summarized using standard goodness-of-fit statistics: root mean squared error (RMSE), normalized RMSE (nRMSE), and Willmott’s d-index.

3. Results and Discussion

3.1. Multi-Sensor Crop Mask and NDVI Eligibility

Across the Texas High Plains, combining CDL crop masks with the NDVI-based eligibility filter produced consistent, NASS-compatible crop areas for both cotton and winter wheat. Figure 1 illustrates this for one cotton and one wheat example. For Gaines County cotton in 2021 (Landsat NDVI), the initial CDL cotton mask mapped 148,103 ha of cotton pixels (Figure 2A). Applying the NDVI eligibility filter to shrink the CDL mask toward the NASS harvested area removed 30,270 ha (≈20%) of low-signal or misclassified pixels, leaving 117,833 ha of NDVI-active cotton (Figure 2B). In this county-year, the iterative eligibility procedure converged to a seasonal NDVI threshold of 0.310, which defines the NDVI-eligible subset shown in Figure 1B and the ineligible pixels in Figure 1C. The resulting NDVI-eligible area closely matches the NASS harvested area for cotton in that county-year.

For Deaf Smith County winter wheat in 2016 (Sentinel-2 NDVI), the initial CDL wheat mask covered 96,576 ha (Figure 2D). Applying the NDVI eligibility filter to shrink the CDL mask toward the NASS harvested area removed 31,280 ha of low-signal or misclassified pixels, leaving 65,296 ha of NDVI-active wheat (Figure 2E). In this county-year, the iterative eligibility procedure converged to an NDVI threshold of 0.291, which defines the NDVI-eligible subset in Figure 2E and the excluded pixels in Figure 2F, and brought the trimmed mask into close agreement with the NASS harvested wheat area within the user-defined tolerance.

Two points are worth emphasizing: (1) The CDL is not used blindly. Without NDVI eligibility, CDL over- and under-estimation relative to NASS can be substantial and spatially structured. The quantile-based NDVI filter removes the noisiest portions of the CDL mask while re-anchoring the total eligible area to NASS harvested area (Section 2.3.1), so the subsequent disaggregation always starts from a NASS-consistent crop footprint rather than raw CDL pixels. In Texas, a significant portion of winter wheat is managed as dual-purpose wheat, where fields are grazed during vegetative stages and may not ultimately be harvested for grain. Approximately 40-50% of planted winter wheat acreage in Texas is either grazed out, terminated early, or abandoned prior to grain harvest, depending on climatic conditions and market incentives (https://varietytesting.tamu.edu/smallgrains/). The increasing availability of satellite time series, cloud computing platforms, and AI-based modeling tools has accelerated the development of intelligent agricultural systems (Li et al., 2025; Wang et al., 2025; Zhang et al., 2025; Yang et al., 2026). The proposed framework complements these predictive approaches by providing statistically consistent spatial yield baselines that can serve as inputs or validation layers for advanced modeling systems.

Present results indicate that this approach is robust across sensors. Figure 2 shows one Landsat-based cotton example and one Sentinel-2-based wheat example, but analogous MODIS-based NDVI masks (not shown) yielded very similar NDVI-eligible areas in the same county-years. This supports the design choice of allowing sensor flexibility (MODIS, Landsat, Sentinel-2) without materially changing the effective crop area used in the yield disaggregation.

A limitation to acknowledge is that in years or counties with persistent cloud cover, complex mixed pixels, or substantial intercropping, the NDVI filter may remove some genuine crop pixels, particularly for coarser MODIS imagery. The union/intersection rules with neighboring years partly mitigate this, but a formal sensitivity analysis of the NDVI quantile thresholds would be a logical extension of this work.

3.2. NDVI-ET Weight Surfaces

For both crops, the combined NDVI-ET weights reproduced expected within-county gradients in crop vigor and seasonal water use within each county. In the Gaines County cotton example for 2021 (Figure 3), the original NDVI exhibits a heterogeneous mosaic of fields with moderate-high greenness (Figure 3A), reflecting variation in crop condition and management, while seasonal ET highlights zones of consistently higher water use, particularly irrigated blocks (Figure 3C). After min-max scaling, the NDVI-based weight field (Figure 3B) enhances contrast among fields with subtle greenness differences, whereas the ET-based weight (Figure 3D) emphasizes fields with persistently higher water use. The raw combined weight (Figure 3E; 0-1 scale) concentrates mass in pixels that are simultaneously green and water-consuming, representing areas with both strong canopy development and sustained transpiration, and the final normalized weight (Figure 3F; ≈5×10⁻⁷-2×10⁻⁶) further redistributes this mass so that the weights sum to one across the NDVI-eligible cotton area.

A similar pattern is evident for winter wheat in Deaf Smith County in 2016 (Figure S1). Original NDVI and ET again show fine-scale spatial heterogeneity across the eligible wheat area (Figure S1A & C), with NDVI emphasizing canopy density and ET capturing differences in seasonal water use. Their respective weight fields (Figure S1B & D) preserve this heterogeneity while compressing values to [0,1]. The combined weight (Figure S1E) tends to favor wheat fields that are simultaneously greener and more water-consuming during the critical spring growth period, over the spring season, and the final normalized weights (Figure S1F; ≈1×10⁻⁶–3×10⁻⁶) identify a subset of fields and partial fields that collectively account for most of the “probability mass” used in subsequent disaggregation.

Because all weights are min-max scaled and then renormalized to sum to one within the NDVI-eligible mask, they can be interpreted probabilistically: a pixel with weight 2×10⁻⁶ contributes roughly twice as much to the county mean yield as a pixel with weight 1×10⁻⁶. The absolute magnitudes are small because the total mass (1) is distributed over tens of thousands of 30-m pixels; spatial allocation depends solely on relative contrasts. This probabilistic interpretation underpins the later β-exponent step, where the weight field can be sharpened (β>1) or flattened (β<1) without changing the county-level NASS constraint.

A limitation that should be acknowledged is that the ET field from OpenET SSEBop is itself model-derived and can smooth or damp local extremes, particularly in heterogeneous irrigated systems. Similar spatial smoothing effects in satellite-based ET products have been noted in previous remote sensing studies, where field-scale irrigation variability or ET estimates can be dampened by model parameterization and spatial averaging (Liu et al., 2026). Where irrigation management varies strongly within a county, ET-based weights may under-represent the true yield contrast among fields, which in turn limits how much spatial variability the disaggregation can express even when β is increased.

Although combining NDVI with ET partially mitigates this issue by integrating canopy vigor and water-use signals, the quality of the resulting weight surface will still depend on the fidelity of both satellite products in capturing within-county variability. The complementary behavior of vegetation indices and evapotranspiration in capturing biomass development and transpiration-driven productivity has been widely demonstrated in intelligent agricultural monitoring systems (Zhai et al., 2023).

3.3. Segmentation and Harvest-Unit Selection

The workflow converts the NDVI-ET weight surface into spatially contiguous harvest units through a segmentation procedure designed to approximate realistic harvest units while avoiding pixel-level noise. Adjacent NDVI-eligible pixels are grouped into contiguous patches, very small fragments are removed, and large patches are optionally subdivided into tiles to ensure that all resulting units fall within a plausible range for the crop and region. This ensures that subsequent disaggregation operates at a scale consistent with operational harvest or management units rather than individual pixels. Figure 4 illustrates the resulting segmentation for cotton in Gaines County in 2021. The output consists of spatially explicit harvest units (Figure 4A), each characterized by its area and associated aggregated weight density (Figure 4B). A comparable segmentation for winter wheat in Deaf Smith County in 2016 is shown in Figure S2. Aggregating fine-scale information into coarser modeling units is a standard practice in integrated environmental and agricultural modeling, where the goal is to retain dominant system behavior while reducing dimensionality and noise (Argent, 2004).

In the Gaines County cotton example for 2021 (Figure 5), the cumulative weight curve rises gradually: the top 1096 of 1,247 units (87.9%) are sufficient to reach the NASS harvested area, and at this point the cumulative weight has already reached 99.6% (Figure 5A). The selected subset covers 114,119 ha compared with a NASS harvested area of 114,121 ha, demonstrating close area anchoring (Figure 5B). This relatively flat cumulative-weight curve indicates that yield contributions are spatially diffused across many units rather than concentrated in a small core, consistent with a mosaic of irrigated fields interspersed with lower-productivity dryland corners. The log-scale histogram of unit weights spans approximately three orders of magnitude, with many low-weight fringe units and a relatively higher number of high-weight units that carry most of the probability mass (Figure 5C).

In contrast, winter wheat production in Deaf Smith County in 2016 is more concentrated: only 474 of 954 units (49.7%) are required to match the NASS harvested area, yet these units already account for 92.1% of the cumulative weight (Figure S3A & B). This steeper cumulative-weight curve indicates that a comparatively small subset of units dominates production, consistent with large, contiguous wheat blocks with broadly similar management and canopy development. The wheat weight histogram is correspondingly broader, with fewer extremely high- or low-weight units (Figure S3C).

Taken together, these diagnostics demonstrate that the harvest-unit selection remains is tightly anchored to official harvested area while preserving spatial heterogeneity embedded in the NDVI-ET weight surface. In cotton, production is more spatially diffused, so most units are needed to reach the NASS area, and nearly all units carry non-negligible weight. In winter wheat, a relatively small share of field units accounts for most of the weight and harvested area, consistent with concentrated high-yield, high-ET irrigated blocks surrounded by many low-weight dryland corners. The cumulative-weight curves also provide a transparent diagnostic for sensitivity analysis: steeper curves (e.g. under higher β exponents) imply that a smaller fraction of units could, in principle, explain most of the county yield, whereas flatter curves indicate that yield is inherently more spread out and less compressible into a few “core” production zones.

3.4. NASS-Anchored Disaggregation to Production Units and Zonal Clusters

By construction, the disaggregation preserves the NASS county mean: the area-weighted mean of all unit or cluster yields closely matches the reported NASS value (692 vs. 703 kg ha^-1) for Gaines County cotton in 2021 (Figure 6C). This area-constrained redistribution approach is conceptually aligned with mass-preserving downscaling methods used in environmental modeling, where fine-scale allocation is required to remain consistent with aggregated statistics (Maestrini and Basso, 2018a; Metwally et al., 2019). The value of the workflow therefore lies not in reproducing the county mean itself, but in reconstructing the spatial distribution of yield consistent with that constraint. Cluster yields in Gaines span approximately 400–1,300 kg ha⁻¹, with a pronounced upper tail (Figure 6B). High-yield clusters are spatially concentrated in contiguous zones in the south-central and western parts of the county, while lower-yield clusters dominate peripheral areas (Figure 6A), consistent with documented contrasts between irrigated cores and marginal dryland lands. Aggregating harvest units into graph-based zonal clusters thus smooths very fine-scale noise but retains dominant spatial gradients, providing production zones that both honor the NASS constraint and preserve within-county heterogeneity.

For Deaf Smith County winter wheat in 2016 (Figure 7), unit-level yields span from low values along the western and southern margins of the county to markedly higher yields in the central and north-eastern blocks (Figure 7A). The resulting histogram is right-skewed, with many clusters around 2–3 Mg ha⁻¹ and a smaller number of high-yield zones approaching 6 Mg ha⁻¹ (Figure 7B). Aggregating harvest units into connectivity-based production clusters smooths out fragmented small units while preserving these broad contrasts: cluster mean yields still range from <2 to >5 Mg ha⁻¹, and the cluster area-weighted mean remains within ~1% of the NASS county mean (Figure 7C). Importantly, these production clusters operate at spatial scales (~ 10 km) compatible with coarser gridded soil and climate datasets commonly used in process-based modeling (e.g., SoilGrids, gridded weather forcing). In this sense, clustering provides a scale bridge between pixel-level remote sensing signals and meso-scale crop or environmental models. This is consistent with previous work showing that aggregating sub-field yield signals to meso-scale zones can retain most of the agronomically relevant variability while greatly reducing model complexity (Schepers et al., 2004; Xiang et al., 2007; Maestrini and Basso, 2018b).

A key control on the contrast of these yield distributions is the exponent β applied to the unit weights. In the baseline configuration shown in Figure 6 and7, disaggregated yields are proportional to $w_u^{\beta }$and then rescaled so that the area-weighted mean matches the NASS yield. When β = 1, yields follow the NDVI-ET weight surface directly; β > 1 amplifies differences and pushes more area toward low and high tails, and β < 1 compresses the distribution toward the county mean. In our sensitivity analysis (Section 2.4), testing β values from 0.8 to 1.8 showed that moderate β (1.2-1.4) can improve agreement with independent variety-trial data in high-yield environments by allowing a small fraction of units or clusters to carry more than their equal share of the county yield, while still honoring the NASS constraint. Excessively large β values, however, generate unrealistically heavy tails and reduced spatial plausibility. Exponent-based transformations to control contrast and heterogeneity are commonly used in spatial allocation and weighting frameworks to regulate distribution tails while preserving aggregate constraints.

Taken together, these results show that the workflow can translate a single NASS county mean into a spatially explicit harvest units and zonal clusters that (i) preserve the official statistics by construction, (ii) reflect realistic spatial patterns driven by NDVI and ET, and (iii) expose tunable yield contrasts through β. This approach therefore provides a reproducible bridge between county statistics, remote-sensing information, and the spatial scales at which process-based models and soil datasets typically operate.

3.5. Evaluation Against Variety Trials

Independent variety trial data provides a demanding test of the disaggregation, because the trials are not located at the centers of the production clusters and, for some sites, coordinates are only approximate (nearest town or station). Under these conditions, the validation question is deliberately strict: within a given search radius around each trial, does the workflow produce any production cluster whose mean yield is close to the observed trial mean?. In addition, the framework does not aim to predict absolute field-level yields, but rather to reconstruct spatial variability under a statistical constraint. Therefore, validation is conducted at the production-zone scale, which is consistent with the spatial resolution of independent observations and the intended spatial decision-support applications.

Across all cotton and winter wheat site-years, performance improved systematically as the search radius increased from 10 to 30 km and as the β-exponent moved from 0.8 toward moderate contrast (β ~ 1.2-1.6), but the gains were incremental rather than dramatic. With a 10 km radius, RMSE ranged from 373 to 389 kg ha^-1 (nRMSE 15.1-15.7%) and Willmott’s d-index from 0.86 to 0.87 for cotton (Figure 8) and from 896 to 967 kg ha⁻¹ (nRMSE 11.7-12.6 %) and Willmott’s d-index from 0.87 to 0.85 across β values (0.8-1.8) for winter wheat (Figure S4). At 20 km, errors dropped to RMSE 282-312 and 713-751 kg ha⁻¹ (nRMSE 11.4-12.6% and 9.3-9.8%, d 0.90-0.93 and 0.90-0.92) for cotton and winter wheat, respectively. Increasing the radius further to 30 km yielded only modest additional improvement, with RMSE of 265-302 kg ha⁻¹ and nRMSE 10.7-12.2% for cotton, and RMSE of 629-677 kg ha⁻¹ and nRMSE ≈ 8.2-8.8% for winter wheat (d ≈ 0.92). Thus, most of the benefit comes from moving beyond a very tight 10 km search window, while 20–30 km radii give broadly similar levels of agreement.

The spatial maps in 2018-2019 for cotton and in 2018 and 2021 for winter wheat complement these summary statistics by showing how cluster yields and trial sites co-vary across Texas (Figure 9). In both cotton and winter wheat, NASS-anchored clusters (dots) form coherent production zones that align with known regional patterns: higher yields in the northern and central High Plains and lower yields toward drier western or southern environments. Variety-trial sites (plus signs) generally sit within or adjacent to clusters of similar magnitude, rather than appearing as isolated high-yield “spikes” in low-yield surroundings. Years with widespread drought or stress (e.g. 2019 in parts of the Panhandle for cotton and 2018 for winter wheat) show spatially consistent yield reductions in both clusters and trials, whereas better years exhibit broader areas of high-yield clusters with trial means that fall within the local cluster range. For example, winter wheat yield exhibited high values in north and east Texas in 2021 in both clusters and trials (Figure 9). These maps therefore demonstrate that the workflow is not only matching trial means in aggregate but is also reproducing plausible spatial gradients across diverse environments.

Although the framework operates under the additional constraint that the area-weighted mean must equal the NASS county statistic, the achieved nRMSE and d-index values fall within the range commonly reported for regional-scale remote-sensing yield estimation at similar spatial resolutions. The β-exponent exerts a clear but non-monotonic influence. At fixed radius, increasing β from 0.8 to ~1.2 improves fit by allowing high-weight clusters (high NDVI and ET) to deviate further above the county mean. Beyond β~1.6, RMSE and nRMSE rise again and declines, indicating that overly aggressive contrast amplification creates too many very low-yield clusters to compensate for a few extreme high-yield ones. In other words, moderate β values (~1.0-1.2) appear to strike the best balance between expressing realistic yield gradients and respecting the NASS constraint. In the absence of independent validation data, a default value of β ≈ 1 provides a robust baseline that preserves the relative structure of spatial variability without introducing excessive contrast. Sensitivity analysis indicates that spatial patterns remain stable across a reasonable range of β values, allowing the parameter to be interpreted as a contrast-adjustment factor rather than a strictly calibrated parameter.

Qualitatively, the scatterplots reveal two consistent behaviors. First, the highest trial yields (>2.0–2.5 Mg ha⁻¹ for cotton and ~6–7 Mg ha⁻¹ for winter wheat) are generally more difficult to reproduce exactly, even under the best β–radius combinations. This compression is expected because the disaggregation is constrained by the county-level NASS mean and cannot generate many extreme high-yield clusters without introducing compensating low-yield areas elsewhere. Second, within the main yield range represented by most observations (~1.0–1.5 Mg ha⁻¹ for cotton and ~2–5 Mg ha⁻¹ for winter wheat), cluster yields were distributed around the 1:1 line with limited apparent bias. This suggests that the NDVI–ET weighting combined with NASS anchoring captured the dominant spatial contrasts among production environments, including irrigated versus dryland systems and broader north–south productivity gradients. Overall, the comparison indicates that the NASS-anchored zonal clusters provide a plausible representation of production-zone-scale yield variability, while very high-yield experimental sites remain difficult to reproduce under a county-mean constraint.

3.6. Synthesis and Limitations

Taken together, the results show that the proposed nayd workflow provides a transparent and reproducible approach for converting aggregated crop statistics into spatially explicit production units and zonal clusters. Rather than predicting yield independently from many environmental covariates, the workflow uses crop masks, multi-sensor NDVI, satellite-derived evapotranspiration, and official yield statistics within a constrained allocation structure. This design allows within-county yield variability to be reconstructed while preserving consistency with reported harvested area and county mean yield. The NDVI-based eligibility step refines the crop extent toward the harvested area reported by NASS, while the combined NDVI–ET weighting captures dominant spatial gradients in crop vigor and water use.

The main contribution of the workflow is that satellite-derived indicators are embedded directly within a statistically consistent spatial allocation process. The resulting production zones are therefore not unconstrained yield predictions, but NASS-consistent representations of crop production patterns. This distinction is important for sustainability assessment and regional decision support, where spatially detailed information is needed but field-level yield observations, proprietary yield-monitor data, or detailed management records are often unavailable. The present implementation deliberately relies on widely available public datasets to maximize reproducibility and scalability across crops, years, and regions.

Validation against independent variety-trial data showed that the resulting production zones reproduced broad site-level yield patterns at approximately 10–20 km scales. The β-exponent provided a transparent mechanism for controlling spatial contrast while preserving the county-level yield constraint, with values near β ≈ 1 offering a stable default configuration when independent calibration data are unavailable.

Importantly, the framework is not restricted to USDA-NASS or CDL products. The same logic can be applied wherever analogous crop masks, aggregated yield statistics, and satellite-derived vegetation or water-use indicators are available. The R-based implementation also allows users to substitute region-specific datasets, adjust growing-season windows, modify eligibility thresholds, and tune clustering scales according to the needs of crop modeling, environmental assessment, or sustainable agricultural planning. Overall, the results demonstrate that combining official agricultural statistics with remote-sensing indicators can generate reproducible and statistically consistent spatial yield layers for applications that require information beyond coarse administrative units.

4. Conclusions

This study developed nayd, a reproducible geospatial workflow for converting aggregated county-level crop statistics into spatially explicit production units and zonal clusters using publicly available crop-mask, NDVI, and evapotranspiration datasets. Unlike conventional remote-sensing yield-mapping approaches that estimate yield independently from satellite predictors, nayd treats official harvested area and county mean yield as fixed statistical constraints. Satellite-derived NDVI and evapotranspiration are used to guide the spatial allocation of those statistics, while NDVI-based crop eligibility, production-unit segmentation, graph-based clustering, and β-controlled allocation provide a transparent structure for reconstructing within-county yield variability without breaking consistency with reported aggregate records.

Application to cotton and winter wheat in Texas showed that the framework can refine CDL crop masks toward reported harvested area, preserve realistic spatial gradients associated with irrigated and dryland systems, and generate production zones compatible with regional crop and soil modeling scales. The β parameter provides an interpretable mechanism for controlling spatial contrast without violating the county-level yield constraint, with values near unity offering a stable default configuration when independent calibration data are unavailable. Evaluation against independent variety-trial data indicated that the resulting NASS-anchored production zones reproduce broad yield patterns at spatial scales relevant for regional analysis and management, rather than attempting to represent exact field-level performance. More broadly, nayd addresses a practical information gap in agricultural data systems: official statistics provide trusted and temporally consistent yield records but are spatially coarse, whereas satellite observations capture spatial heterogeneity but are not inherently linked to reported production totals. The workflow connects these two information streams in a form that is directly usable for regional crop modeling, environmental assessment, geospatial analysis, and agricultural decision support. It can also be adapted to other crops and regions where analogous datasets are available.