Grain yield maps must accurately display general yield patterns as well as details of local yield variation. Different geostatistical procedures for creating interpolated yield maps by integrating yield data with remotely sensed vegetation indices (VI) were evaluated. Yield monitor data and a multispectral satellite image at 4-m spatial resolution were collected in an irrigated maize (Zea mays L.) field and a rainfed soybean [Glycine max (L.) Merr.] field. Ordinary kriging (OK), cokriging (CK), simple kriging with varying local means (SKLM), and kriging with external drift (KED) were compared using grain yield as primary variable and three different VI as secondary variables. At both sites, SKLM performed best in terms of the precision of grain yield maps and maps that depicted true yield patterns. Utilizing the most suitable vegetation index at each site and SKLM as interpolation method, the root mean squared error of yield predictions was increased by nearly 20% over OK. Methods such as KED or CK resulted in only small improvement of yield maps over those obtained with OK. Utilizing the satellite image decreased errors associated with yield monitor operations and allowed better prediction in areas where no reliable yield measurements were available. Due to the robust nature of modeling the relationship between primary and one or more secondary variables, the SKLM method has considerable potential for use in commercial precision-farming software. Future efforts on improving yield mapping should concentrate on obtaining improved yield monitor data and imagery and developing yield-sensitive VI for high-yielding crops.