National Mapping of New Zealand Pasture Productivity Using Temporal Sentinel-2 Data

: A national map of pasture productivity, in terms of mass of dry matter yield per unit area and time, enables evaluation of regional and local land-use suitability. Difﬁculty in measuring this quantity at scale directed this research, which utilises four years of Sentinel-2 satellite imagery and collected pasture yield measurements to develop a model of pasture productivity. The model uses a Normalised Difference Vegetation Index (NDVI), with spatio-temporal segmentation and averaging, to estimate mean annual pasture productivity across all of New Zealand’s grasslands with a standard error of prediction of 2.2 t/ha/y. Regional aggregates of pasture yield demonstrate expected spatial variations. The pasture productivity map may be used to classify grasslands objectively into stratiﬁed levels of production on a national scale. Due to its ability to highlight areas of land use intensiﬁcation suitability, the national map of pasture productivity is of value to landowners, land users, and environmental scientists.


Introduction
New Zealand agricultural land, which historically has primarily been used for sheep grazing, has seen increased utilisation for dairy in recent decades. Indeed, the land area of dairy farming increased by 33% between 1999 and 2019 [1,2]. Land use intensification in dairy farming areas has also occurred, with an increase from 1.9 to 2.8 cows per hectare between 1991 and 2019 [2,3]. There have been regional patterns to New Zealand dairy conversions, typically based on the suitability of land for cattle grazing. For instance, Southland has seen an increase in dairy farming land from 65,000 ha in 1999 to 221,000 ha in 2019 [1,2].
One metric of interest in the prediction of land suitability for high-intensity farming is pasture productivity, or pasture yield, expressed in units of tonnes per hectare per year (t/ha/y). Pasture productivity is the mean annual mass of dry matter obtained from a crop per unit area, and thus represents the quantity of feed available for grazing in that location. Understanding the spatial variations of pasture productivity across the New Zealand landscape can provide information on the suitability of land for intensification.
Data for a national map of pasture productivity have proved difficult to obtain [4]. Direct measurement of pasture biomass through mechanical clipping of quadrats has been a typical approach for gathering data at plot scale [5], which may be extended to field scale with capacitance probes [6]. While plot-based measurements are accurate, repetitive measurements over large areas are expensive, so remote sensing is required to scale up to regional or national scale [7].
Reinermann et al. [7] recently reviewed 253 research articles on remote sensing of grassland production. They found that utilisation of satellite information had high value for large or remote areas, and that vegetation indices were commonly used in empirical relationships with biomass measurements, with NDVI [8] used in 62% of studies, EVI [9] in 15%, and SAVI [10] in 9%. The R 2 values of biomass models ranged from 0.4 to 0.97, with the highest R 2 obtained when temporal or spatial heterogeneity was restricted (for example, Edirisinghe et al. [4] obtained an R 2 of 0.74 for fields on one farm in October). Long-term biomass estimates have been made using satellite imagery, but these estimates have so far only been applied at a field scale [11], or for single seasons without the ability to represent temporally averaged pasture yield [12].
There is therefore a gap in contemporary research for national mapping of annualised pasture productivity, a metric that would enable the geographic variations of land suitability for intensification to be objectively quantified. Here, we present a method to estimate annualised pasture yield across all New Zealand's grasslands. It uses a long-term Sentinel-2 imagery archive and selected analyses of dry matter yield. Pasture yield is regressed against a temporal median of NDVI in cloud-free pixels. The regression model is applied to a vector layer of land parcel segments, generated using an automated spatio-temporal approach, to estimate spatial variations in pasture productivity across all New Zealand's grasslands. This research provides a proof of concept of a method of using satellite imagery to determine annual pasture productivity across New Zealand.

Development of Pasture Productivity Model from Sentinel-2 Imagery
Pasture yield data were sought in published literature for the timeframe for which Sentinel-2 imagery is available (September 2015-January 2020). However, limited data with accurate georeferencing were available. The locations for which we have data are listed in Table 1. The data obtained from these sources were mean annual yields of tonnes per hectare (t/ha/y). While georeferencing was available as point data, the point-based figures were assumed to be representative of the fields they were within. It was also assumed that management protocols had not changed since yield recordings were obtained. Temporal stacks of annual satellite imagery mosaics were created for each site, extracted from an analysis-ready Sentinel-2 dataset. The images had been top-of-atmosphere corrected using the 6S code [22] and projected to the New Zealand Transverse Mercator (NZTM) map grid. Additionally, the images had been topographically flattened [23], using a digital elevation model (DEM) derived from 20 m contours to normalize brightness variations in steep terrain.
Visual interpretation of the image mosaics was used to generate a polygon just inside the field boundary containing each site, in order to minimise edge effects on polygon means. The polygons were created to take advantage of the homogeneity of individually managed fields and reduce the noise associated with single pixels [24]. Examples of these polygons and the associated point locations are shown in Figure 1. Temporal stacks of annual satellite imagery mosaics were created for each tracted from an analysis-ready Sentinel-2 dataset. The images had been top-of-atm corrected using the 6S code [22] and projected to the New Zealand Transverse (NZTM) map grid. Additionally, the images had been topographically flattened ing a digital elevation model (DEM) derived from 20 m contours to normalize b variations in steep terrain.
Visual interpretation of the image mosaics was used to generate a polygon ju the field boundary containing each site, in order to minimise edge effects on means. The polygons were created to take advantage of the homogeneity of ind managed fields and reduce the noise associated with single pixels [24]. Example polygons and the associated point locations are shown in Figure 1.
We assume that there exists a correlation between the annual yield, and the response of the site. Namely, the mass of dry matter should be related to the lea vegetation, measured as that location's NDVI. For the interpretation of NDVI fro nel-2 imagery in this paper, the relationship is defined as:  We assume that there exists a correlation between the annual yield, and the spectral response of the site. Namely, the mass of dry matter should be related to the leafiness of vegetation, measured as that location's NDVI. For the interpretation of NDVI from Sentinel-2 imagery in this paper, the relationship is defined as: Here, NIR is the top-of-atmosphere corrected intensity of near-infrared reflectance in band 8 of Sentinel-2, and Red is the intensity of reflectance in band 4. Due to the difference Remote Sens. 2021, 13, 1481 4 of 16 in scale of these bands, a factor of four was chosen to scale red reflectance to allow changes in red reflectance to have a similar weight to NIR changes in the NDVI range of interest.
Behaviour indicative of high-producing pasture is quick to return to high levels of vegetation after grazing. Given that continuous imagery of the sites is unavailable, the spatial mean of median NDVI values of each pixel was thought to be an indicator of pasture quality. For each pixel in each polygon, the images that contained valid data (i.e., omitting cloud, cloud shadow, and snow) were sorted by NDVI, and the median taken. The spatial mean of these values within a polygon was assigned as the polygon's median NDVI.
These median NDVIs were plotted against yield and a linear model fitted to the data. A linear model was chosen because the domain of the input data (moderate-to-high producing grassland) was most likely to be linearly related to NDVI. Additionally, use of a more complex model with the limited available data would have created an overfitted model. The generated equation was used, in combination with a national median NDVI image, to create a national map of mean annual pasture productivity.

Generation of a Median NDVI Image of New Zealand
To apply the generated equation to a vector map of pasture fields, a national raster map of NDVI was required. This process required a combined multi-temporal and spatial approach to analysis-ready Sentinel-2 satellite imagery. For each topographically flattened image in the archived sequence, the NDVI was computed using (1). The median NDVI was then computed for each pixel, with dates being omitted from sequences if cloud, cloud shadow, or snow were present on that pixel.
Computation of pixelwise median NDVI for the study sites was straight-forward as they had been visually verified to be all pasture, so simply taking the median NDVI from the cloud-free dates for each pixel was suitable. Extending this method to a national median NDVI map required accounting for seasonal crops and crop rotations that would arbitrarily reduce NDVI in particular images. To produce a median NDVI image which reflected the NDVI of planted crops, the process depicted in Figure 2 was followed. Dates with negative NDVI were assumed to be unvegetated and were therefore ignored in median calculation. The exclusion of these dates is deliberate, because dates with negative NDVI are not representative of the typical pastural vigour of that location, so should not be used to draw conclusions about pasture yield of that location. Here, NIR is the top-of-atmosphere corrected intensity of near-infrared reflectance in band 8 of Sentinel-2, and Red is the intensity of reflectance in band 4. Due to the difference in scale of these bands, a factor of four was chosen to scale red reflectance to allow changes in red reflectance to have a similar weight to NIR changes in the NDVI range of interest.
Behaviour indicative of high-producing pasture is quick to return to high levels of vegetation after grazing. Given that continuous imagery of the sites is unavailable, the spatial mean of median NDVI values of each pixel was thought to be an indicator of pasture quality. For each pixel in each polygon, the images that contained valid data (i.e., omitting cloud, cloud shadow, and snow) were sorted by NDVI, and the median taken. The spatial mean of these values within a polygon was assigned as the polygon's median NDVI.
These median NDVIs were plotted against yield and a linear model fitted to the data. A linear model was chosen because the domain of the input data (moderate-to-high producing grassland) was most likely to be linearly related to NDVI. Additionally, use of a more complex model with the limited available data would have created an overfitted model. The generated equation was used, in combination with a national median NDVI image, to create a national map of mean annual pasture productivity.

Generation of a Median NDVI Image of New Zealand
To apply the generated equation to a vector map of pasture fields, a national raster map of NDVI was required. This process required a combined multi-temporal and spatial approach to analysis-ready Sentinel-2 satellite imagery. For each topographically flattened image in the archived sequence, the NDVI was computed using (1). The median NDVI was then computed for each pixel, with dates being omitted from sequences if cloud, cloud shadow, or snow were present on that pixel.
Computation of pixelwise median NDVI for the study sites was straight-forward as they had been visually verified to be all pasture, so simply taking the median NDVI from the cloud-free dates for each pixel was suitable. Extending this method to a national median NDVI map required accounting for seasonal crops and crop rotations that would arbitrarily reduce NDVI in particular images. To produce a median NDVI image which reflected the NDVI of planted crops, the process depicted in Figure 2 was followed. Dates with negative NDVI were assumed to be unvegetated and were therefore ignored in median calculation. The exclusion of these dates is deliberate, because dates with negative NDVI are not representative of the typical pastural vigour of that location, so should not be used to draw conclusions about pasture yield of that location.  An additional requirement for producing a median NDVI image across New Zealand was the automatic omittance of images with cloud, cloud shadow, and snow for each pixel. To achieve the set of unobstructed views of a pixel as shown in Figure 2, an automated process was utilised to identify scenes for which the pixel was obscured by these phenomena. These conditions were found using single image cloud detection with FMASK [25] feeding into a multitemporal algorithm TMASK [26]. Pixels with values significantly outside the expected range of reflectances were thus ignored [27].

Generation of a National Vector Map of Average Median NDVI
Taking a raster map of median NDVI as an input for a pasture productivity model would result in a noisy output. Pasture is managed at a field level, so mapping productivity to smaller units is not a suitable interpretation of the available data. Instead, as was performed for the study sites, the entire field can be taken as a single unit. This is a more robust approach than analysing data at an individual pixel level [24]. We used the multitemporal method of North et al. [28] to segment the raster imagery into field boundaries using a sequence of imagery from September 2017 to November 2018. It is assumed that these boundaries are suitable divisions of uniform management for the duration of the study. Examples of the output of this process in areas of pasture are shown in Figure 3. An additional requirement for producing a median NDVI image across New Zealand was the automatic omittance of images with cloud, cloud shadow, and snow for each pixel. To achieve the set of unobstructed views of a pixel as shown in Figure 2, an automated process was utilised to identify scenes for which the pixel was obscured by these phenomena. These conditions were found using single image cloud detection with FMASK [25] feeding into a multitemporal algorithm TMASK [26]. Pixels with values significantly outside the expected range of reflectances were thus ignored [27].

Generation of a National Vector Map of Average Median NDVI
Taking a raster map of median NDVI as an input for a pasture productivity model would result in a noisy output. Pasture is managed at a field level, so mapping productivity to smaller units is not a suitable interpretation of the available data. Instead, as was performed for the study sites, the entire field can be taken as a single unit. This is a more robust approach than analysing data at an individual pixel level [24]. We used the multitemporal method of North et al. [28] to segment the raster imagery into field boundaries using a sequence of imagery from September 2017 to November 2018. It is assumed that these boundaries are suitable divisions of uniform management for the duration of the study. Examples of the output of this process in areas of pasture are shown in Figure 3.  The polygons created using this method covered the entirety of New Zealand's agricultural land. To ensure that only pasture crops were assessed for their estimated yield, this vector map was clipped to the high and low producing pasture classes in the Land Cover Database (LCDB) [29]. High producing pasture comprises exotic grassland characterised by a spectral signature indicating good vegetation vigour. It is typically subject to pasture renewal every 5-10 years and consists of species such as clovers, ryegrass, and cocksfoot. Low producing pasture consists of a mixture of exotic (including browntop and sweet vernal) and indigenous (including hard, blue, and silver tussock) grasslands that typically brown off during summer months. The intersection of these two LCDB pasture classes and field polygons with area greater than 0.1 ha was obtained to produce a set of field polygons which could have their pasture yield estimated. Small field polygons were excluded because they were unlikely to represent pasture.
For each of the pasture fields, the pixels with centroids that fell within the polygon boundaries were collected. The mean value of temporal median NDVI (negative NDVI values were excluded) of these pixels was assigned to each polygon. This value was then used as the input to the model equation allowing a pasture yield in t/ha/y to be assigned. Polygons with mean temporal median NDVI outside the domain of the data used to formulate the model were excluded from this calculation and set to have zero pasture yield.

Results
Using the data in Table 1, the model relating median NDVI to mean annual pasture productivity was found to be: The correlation coefficient for this model was 0.70. A plot depicting the pasture yield data with the fitted model relationship is shown in Figure 4. The polygons created using this method covered the entirety of New Zealand's agricultural land. To ensure that only pasture crops were assessed for their estimated yield, this vector map was clipped to the high and low producing pasture classes in the Land Cover Database (LCDB) [29]. High producing pasture comprises exotic grassland characterised by a spectral signature indicating good vegetation vigour. It is typically subject to pasture renewal every 5-10 years and consists of species such as clovers, ryegrass, and cocksfoot. Low producing pasture consists of a mixture of exotic (including browntop and sweet vernal) and indigenous (including hard, blue, and silver tussock) grasslands that typically brown off during summer months. The intersection of these two LCDB pasture classes and field polygons with area greater than 0.1 ha was obtained to produce a set of field polygons which could have their pasture yield estimated. Small field polygons were excluded because they were unlikely to represent pasture.
For each of the pasture fields, the pixels with centroids that fell within the polygon boundaries were collected. The mean value of temporal median NDVI (negative NDVI values were excluded) of these pixels was assigned to each polygon. This value was then used as the input to the model equation allowing a pasture yield in t/ha/y to be assigned. Polygons with mean temporal median NDVI outside the domain of the data used to formulate the model were excluded from this calculation and set to have zero pasture yield.

Results
Using the data in Table 1, the model relating median NDVI to mean annual pasture productivity was found to be: The correlation coefficient for this model was 0.70. A plot depicting the pasture yield data with the fitted model relationship is shown in Figure 4. To generate a map of pasture productivity from this model, a national map of median NDVI was required. Figure 5 shows a map of New Zealand with the pixelwise median of all cloud-free images in the date range November 2015-February 2020. Figure 6 shows a To generate a map of pasture productivity from this model, a national map of median NDVI was required. Figure 5 shows a map of New Zealand with the pixelwise median of all cloud-free images in the date range November 2015-February 2020. Figure 6 shows types, rather than climatic conditions, is visible. Small polygons not likely to be representative of fields, in addition to features such as roads and waterways, have been omitted from the map.
The production of a national map of pasture productivity enables evaluation of regional differences in pasture yield. These data are presented in Table 2. Historically, Waikato and Taranaki's agricultural land has been heavily utilised by New Zealand's dairy industry, so it is thus unsurprising that those regions have the highest mean pasture production (9.8 t/ha/y and 9.9 t/ha/y, respectively) among North Island regions. Additionally, notable are the high average pasture yield figures in the West Coast and Southland regions (10.7 t/ha/y each); these are among New Zealand's wettest and coldest regions. Figure 5. A national map of median NDVI, with regional boundaries drawn in black. Figure 5. A national map of median NDVI, with regional boundaries drawn in black.
Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 16 Figure 6. A national map of median NDVI, omitting dates for each pixel where the NDVI is negative. Regional boundaries are shown in black. Figure 6. A national map of median NDVI, omitting dates for each pixel where the NDVI is negative. Regional boundaries are shown in black. Figure 6 was used in conjunction with the developed model to create a national pasture productivity map of New Zealand. Figure 7 shows the yield in areas planted with pasture (according to LCDB) down to NDVI = 0.2 (lower limit of predicted yield is therefore 1.7 t/ha/y). To complement the national level maps, regional maps have also been produced to demonstrate spatial variations in pasture productivity. For example, Figure 8 depicts the Otago region's pasture yield, with below-threshold production grasslands in black.    This image demonstrates expected localised differences in yield. There is notably higher production in wetter areas, and lower production in the drier areas, where more marginal grasslands are found. Similar patterns are observed in Taranaki, shown in Figure 9. There are expected patterns of high-and low-producing grassland visible. Figure 8. Otago pasture yield, with black areas having lower productivity than can by this model. Figure 9. Taranaki pasture yield, with black areas having lower productivity than c by this model. Figure 9. Taranaki pasture yield, with black areas having lower productivity than can be predicted by this model. Figure 10 depicts pasture yield in mid-Canterbury. At this scale, the variability of pasture yield due to management strategies and other localised differences such as soil types, rather than climatic conditions, is visible. Small polygons not likely to be representative of fields, in addition to features such as roads and waterways, have been omitted from the map.
The production of a national map of pasture productivity enables evaluation of regional differences in pasture yield. These data are presented in Table 2. Historically, Waikato and Taranaki's agricultural land has been heavily utilised by New Zealand's dairy industry, so it is thus unsurprising that those regions have the highest mean pasture production (9.8 t/ha/y and 9.9 t/ha/y, respectively) among North Island regions. Additionally, notable are the high average pasture yield figures in the West Coast and Southland regions (10.7 t/ha/y each); these are among New Zealand's wettest and coldest regions.

Accuracy of Results
The paucity of field data available for model development meant that it was infeasible to remove a random selection of data points for accuracy evaluation without adversely affecting the quality of the model. We thus assumed that the data in Table 1 were randomly spatially sampled to evaluate the accuracy of the map. This assumption enabled the calculation of the standard error of prediction [30] for the model to be calculated while maintaining the integrity of the model. By this means, the standard error was calculated to be 2.2 t/ha/y at 70% confidence level.
The data collated in Table 1 were suited to a linear model and had a range of pasture yields from 4.8 t/ha/y to 17.2 t/ha/y. Median NDVI in the corresponding range of the input data could be used in the model to estimate pasture yield. Outside the range of measured pasture yields, this model is less applicable. For instance, the linear trendline indicated that pasture yield would be nil if the mean median NDVI dropped below 0.157; however visual inspection of sites with values below this showed pasture which must have some positive yield. Therefore, it is likely that pasture yield with respect to median NDVI is asymptotic as it approaches 0 t/ha/y. More field data from sites with very low pasture yield are needed to extend the model. Seasonal effects on pasture productivity could also be assessed if seasonlized field data was available.
The data collated in Table 1 were provided with varying methods of error assessment, and some data were missing error assessment entirely. While efforts were made to obtain accurate data from sources with known quality, it was not possible to standardize the error of the input data. Future studies may improve on the method presented in this paper by biasing the model regression with more accurately measured data.

Computational Complexities
Generating a national map of pasture productivity was computationally expensive. The volume of data and duration of processing necessitated the use of New Zealand's National eScience Infrastructure (NeSI), so that we could leverage the parallelisable elements of the experimental procedure. The step of computing pixelwise median NDVI required the sorting of up to 120 values for each 10 m pixel and was computed separately on each of the five satellite passes over New Zealand. The agricultural land area then needed to be segmented using the temporal method detailed. This was achieved by separating the landscape according to New Zealand's 16 regions [31] and running the segmentation algorithm with as single CPU devoted to each region. This process took 24 h for the largest region (Canterbury). Separation of the national land parcel polygon file into regional files circumvented file size limit issues related to the shapefile format. All these considerations were important due to the volume of the data processing required.

Implications of LCDB
The national pasture productivity map was evaluated over polygons from LCDB with 2012 and 2018 class associations 40 (High Producing Exotic Grassland) and 41 (Low Producing Grassland). The union of these polygons was taken as New Zealand's grassland area, without consideration taken to the previously determined production category. The differences between these two classes, according to the class definition documentation [29], is according to a combination of species and the visibility (or lack thereof) of "lush growth due to inherently high soil fertility or high annual rainfall". Furthermore, the analysis of this growth can be reflective of the dates of image acquisition. These processes are subjective and potentially prone to bias, so this national map is a means to provide an objective measure of the relative productivity of pasture areas. Figure 11 depicts the distribution of high producing and low producing grasslands according to LCDB. Overlaid onto this map are areas that are (a) high producing as classed in LCDB, but have mean median NDVI below 0.2, and (b) low producing according to LCDB, but have mean median NDVI above 0.2. This threshold, which, using the model defined in "Results", corresponds to a yield of 1.7 t/ha/y, was selected as a best guess according to expert remote sensing knowledge. At a national scale, this threshold appears to reflect the delineations set by LCDB mappers quite well. However, disagreements in classification are more clearly visible at finer scales. An example of this is shown in Figure 12.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 16 prioritise (or otherwise) particular land for dairy intensification while minimising negative environmental impacts. Figure 11. National map of high producing (pink) and low producing (blue) grasslands. Darker shades indicate areas which have been misclassified according to the median NDVI method presented in this paper.     To assess the accuracy of the LCDB classifications against this objective metric, a balanced accuracy approach was used. The two classes are uneven in terms of area, with 7.6 million hectares of land classified as high producing by LCDB, as opposed to 1.4 million hectares of low producing grassland. This class disparity meant that the average accuracy needed to be assessed for each class. With a threshold of 0.2, 81.4% of the balanced area of grassland is classified in agreement with the objective method. Further investigation of the remaining 18.6% of grassland area is required to determine whether an objective definition would improve the classification between high and low producing grassland.

Further Implications for Science
The method of temporal analysis presented in this paper has further applications in identification of land cover types. The median NDVI approach allows for analysis of land cover over long temporal sequences, with seasonal changes in spectral response appropriately smoothed out. Furthermore, the identification of non-vegetated images in pixelwise sequences allows for cropping cycles to be investigated. As explained in the methods section, when estimating pasture yield, non-vegetated images should be ignored. However, the frequency of such images in a sequence can be used to identify cropping land. Figure 13 shows an example of this. The ratio of the number of vegetated observations of a pixel to the total number of observations of that pixel was calculated nationally for the time series of Sentinel-2 images. A three-band RGB stretch of median NDVI of vegetated observations, median NDVI of all observations, and vegetated ratio is shown in the left image. This image highlights areas where short-rotation crops are planted (showing as red in the image), because these areas have varying NDVI levels over time, and in some instances have negative NDVI values which are ignored for the calculation of vegetated median NDVI. These areas are in contrast to areas showing as blue; these have mostly identical NDVI time series and therefore have similar red and green band values. Distinguishing these areas of cropland as opposed to permanent pasture using single images or mosaics can be comparatively difficult, depending on the date of image acquisition and planted crop. Thus, this objective method can be used to identify areas in the LCDB where, although land is marked as containing pasture, the multitemporal spectral response of those areas' pixels indicates that the land cover is on a rotation.  Another application of this approach can be seen in Figure 14. The LCDB is derived from national mosaics of summer imagery. Distinguishing areas populated with indigenous evergreen species from deciduous species in these summer mosaics is difficult because their spectral responses in summer are so similar. This similarity can be seen in the darker areas of forest in the right-hand side image. The NDVI median stretch shown on the left highlights the deciduous areas in red, and evergreen areas in white. These groupings of species constitute distinct cover classes in LCDB, so this temporal approach can assist in differentiating evergreen from deciduous cover classes.

Conclusions
A national map of mean annual pasture productivity has been produced from a timeseries of Sentinel-2 satellite imagery. The linear model used to produce this map relates median NDVI to mean annual pasture yield with a standard error of prediction of 2.2 t/ha/y. The model was applied to field-level segments, enabling regional and local trends to be explored. The map will be used in the future to objectively classify regions of pasture The national map of pasture productivity provides information about the suitability of land for specific agricultural purposes. Areas of land with greater ability to yield dry matter for stock feed are more suitable for high intensity dairy farming than areas with lower pasture yield. Modelling of this quantity, without the ability to carry out the expensive and time-consuming process of directly measuring yield [5], would typically require prior knowledge of intraregional variation of soil types in combination with geographic slope and aspect information. The approach presented in this paper circumvents the need to obtain these data by evaluating vegetation response directly. As such, the combination of this map and knowledge of local land management strategies can help stakeholders prioritise (or otherwise) particular land for dairy intensification while minimising negative environmental impacts.

Conclusions
A national map of mean annual pasture productivity has been produced from a time-series of Sentinel-2 satellite imagery. The linear model used to produce this map relates median NDVI to mean annual pasture yield with a standard error of prediction of 2.2 t/ha/y. The model was applied to field-level segments, enabling regional and local trends to be explored. The map will be used in the future to objectively classify regions of pasture by their productivity. This temporal method of NDVI analysis can also be used to differentiate between permanent pasture and short-rotation cropland, and also between evergreen and deciduous cover classes.

Conflicts of Interest:
The authors declare no conflict of interest.