Gravitational wave (GW) observatories are discovering binary neutron star mergers (BNSMs), and in at least one event we were able to track it down in multiple wavelengths of light, which allowed us to identify the host galaxy. Using a catalogue of local galaxies with inferred star formation histories and adopting a BNSM delay time distribution (DTD) model, we investigate the dependence of BNSM rate on an array of galaxy properties. Compared to the intrinsic property distribution of galaxies, that of BNSM host galaxies is skewed towards galaxies with redder colour, lower specific star formation rate, higher luminosity, and higher stellar mass, reflecting the tendency of higher BNSM rates in more massive galaxies. We introduce a formalism to efficiently make forecast on using host galaxy properties to constrain DTD models. We find comparable constraints from the dependence of BNSM occurrence distribution on galaxy colour, specific star formation rate, and stellar mass, all better than those from dependence on r-band luminosity. The tightest constraints come from using individual star formation histories of host galaxies, which reduces the uncertainties on DTD parameters by a factor of three or more. Substantially different DTD models can be differentiated with about 10 BNSM detections. To constrain DTD parameters at 10 per cent precision level requires about one hundred detections, achievable with GW observations on a decade time-scale.