I have tried calculating and plotting these with HiEst, Introgress etc but with no luck. Is there a more modern method to do this, using ddRADseq data?
Yep - that's right Ricardo. I ended up using the built-in method for estimating interspecific heterozygosity and hybrid index in bgc, but I'm not sure it's as good as other methods.
Introgress has a triangle.plot() function. The hi.index and int.het could come from any source (i.e. they don't necessarily need to be calculated in Introgress).
Also of note, a colleague and I have some R code for producing other plots from bgc outputs, since it sounds like you've run that already: https://github.com/btmartin721/ClinePlotR
I've also got some scripts for creating inputs for Introgress and bgc that you might find helpful: https://github.com/tkchafin/scripts (see phylip2bgc.pl and phylip2introgress.pl). The intended input for both is a phylip file of concatenated SNPs (e.g. from ddRAD data processed in ipyrad). The ClinePlotR package also has some functions for file conversions.