I intend to fit spatial covariance structure models such as spherical, Gaussian, exponential, power, linear, and linear log to a yield monitor data in order to identify the structure that best describes the spatial patterns.

Yield monitor data is a kind of spatial data obtained from combine harvester equipped with yield monitor device linked with GPS. Here, the yield is the attribute component georeferenced with longitude and latitude in decimal degree.

Firstly, I start with exploratory spatial data analysis by plotting the spatial locations of sampling point having converted the spherical coordinates into easting and northing projected coordinates. Being a massive dataset with 12,438 observations, I find it difficult to visualize whether there is a presence of surface trends or not before proceeding with the actually modeling.

Could anyone guide me on how to know if there is trend in this situation with large dataset? If there is a trend, I have to model the residual obtained from regressing the yield as a function of spatial coordinates but now I never know whether there is trend or not.

Is it proper to possibly divide the data into certain number of subregion and fit the spatial structure to each region by having small number of observations? Please help me out. See the attached exploratory analysis.