We introduce a new effective strategy to assign group and cluster membership probabilities $P_{mem}$ to galaxies using photometric redshift information. Large dynamical ranges both in halo mass and cosmic time are considered. The method takes the magnitude distribution of both cluster and field galaxies as well as the radial distribution of galaxies in clusters into account using a non-parametric formalism and relies on Bayesian inference to take photometric redshift uncertainties into account. We successfully test the method against 1,208 galaxy clusters within redshifts $z=0.05-2.58$ and masses $10^{13.29-14.80}~M_odot$ drawn from wide field simulated galaxy mock catalogs developed for the Euclid mission. Median purity $(55^{+17}_{-15})%$ and completeness $(95^{+5}_{-10})%$ are reached for galaxies brighter than 0.25$L_ast$ within $r_{200}$ of each simulated halo and for a statistical photometric redshift accuracy $sigma((z_s-z_p)/(1+z_s))=0.03$. The mean values $p=56%$ and $c=93%$ have sub-percent uncertainties. Accurate photometric redshifts ($sigma((z_s-z_p)/(1+z_s))lesssim0.05$) and robust estimates for the cluster redshift and the center coordinates are required. The method is applied to derive accurate richness estimates. A statistical comparison between the true ($N_{rm true}$) vs estimated richness ($lambda=sum P_{mem}$) yields on average to unbiased results, $Log(lambda/N_{rm true})=-0.0051pm0.15$. The scatter around the mean of the logarithmic difference between $lambda$ and the halo mass is 0.10~dex, for massive halos $gtrsim10^{14.5}~M_odot$. Our estimates could be useful to calibrate independent cluster mass estimates from weak lensing, SZ, and X-ray studies. Our method can be applied to any list of galaxy clusters or groups in both present and forthcoming surveys such as SDSS, CFHTLS, DES, LSST, and Euclid.