In this article, we present an approach based on generalized additive models (GAMs) to predict species' distributions and abundance in Florida estuaries with habitat suitability modeling. Environmental data gathered by fisheries-independent monitoring in Tampa Bay from 1998 to 2008 were interpolated to create seasonal habitat maps for temperature, salinity, and dissolved oxygen and annual maps for depth and bottom type. We used delta-GAM models assuming either zero-adjusted gamma or beta-inflated-at-zero distributions to predict catch per unit effort (CPUE) from five habitat variables plus gear type for each estuarine species by life stage and season. Bottom type and gear type were treated as categorical predictors with reference parameterization. Three spline-fitting procedures (the penalized B-spline, cubic smoothing spline, and restricted cubic spline) were applied to the continuous predictors. Two additive, linear submodels on the log and logistic scales were used to predict CPUEs >0 and CPUEs = 0, respectively, across environmental gradients. The best overall model among those estimated was identified based on the lowest Akaike information criterion. A stepwise routine was used to omit continuous predictors that had little predictive power. The model developed was then applied to interpolated habitat data to predict CPUEs across the estuary using GIS. The statistical models, coupled with the use of GIS, provide a method for predicting spatial distributions and population numbers of estuarine species' life stages. An example is presented for juvenile pink shrimp Farfantepenaeus duorarum during the summer in Tampa Bay, Florida.