We consider the problem of estimation of a bivariate density function with support $\Re\times[0,\infty)$, where a classical bivariate kernel estimator causes boundary bias due to the non-negative variable. To overcome this problem, we propose four kernel density estimators and compare their performances in terms of the mean integrated squared error. Simulation study shows that the estimator based on the proposed normal-gamma kernel performs best. Two astronomical data sets are used to demonstrate the applicability of this estimator.