Your approach seems reasonable, especially your choice to stratify your sampling. This will make it more efficient provided you can easily delineate the different industrial zones.
I don't have a book to recommend you, but you could model your uncertainty using the Poisson distribution, with the $\lambda =$ No. of Toxic Waste Sites per Square Kilometer.
You could carry out your sampling program as you described and then find the maximum likelihood estimator for $\lambda_{Ai}$ where A is the area of a sampling sector in zone $i$. In particular, you would maximize the following formula wrt $\lambda_{Ai}$ where $N_i$= number of sectors sampled from zone $i$:
$\max\limits_{\lambda_{Ai}} \prod\limits_{j=1}^{N_i} \frac{e^{-\lambda_{Ai}}\lambda_{Ai}^{n_{ij}}}{n_{ij}!}$ where $n_{ij}$ is the number of toxic sites in sector $j$ of zone $i$. The value of $\lambda_{Ai}$ that maximizes the product is $\lambda_{Ai}^* = \frac{1}{N_i}\sum\limits_{j=1}^{N_i}{n_{ij}}$
You will get one estimate per zone, $\lambda_{Ai}^*$, which you can interpret as the frequency of toxic waste sites within a region of area $A_i$. Your uncertainty for the total number of sites in Zone $i$ with total area $A_{Ti}$can be modeled using your estimated $\lambda_{Ai}^*$ in the Poisson distribution: $Poisson(\lambda_{Ai}^*\frac{A_{Ti}}{A_i})$.
To get a country-wide estiamte, you would need to combine the $\lambda_{Ai}^*$ into another Poisson distribution: Total No. of Sites ~ $Poisson(\sum\limits_{i=1}^{N_{zones}}\lambda_{Ai}^*\frac{A_{Ti}}{A_i})$.
Refinements
The above should get you a decent estimate. However, if your country is small enough that your sample will cover an appreciable portion of the total land area or area within a zone, then you should reduce the total area for each zone by the sampled area in the above formula, so you are modleing the uncertainty on the remaining area (which is actually more accurate in both cases), then you add this uncertainty to your actual counts in the areas you've sampled.
Also, you will notice that you're using a point estimate of $\lambda$. There is some uncertianty in the actual value of this quantity, but including it requires using extended likelihood for predicting a Poisson variable. The formula is pretty simple, if Y is the total number of sites in zone $i$, then the likelihood function for Y is:
$L(Y_i) = e^{-(N+1)\hat\theta(Y_i)}\frac{\hat\theta(Y_i)^{Y_i+\sum\limits_{j=1}^{N_i}{n_{ij}}}}{Y_i!}$ Where $\hat\theta(Y_i) = \frac{A_{Ti}}{A_i}(Y_i + \frac{\sum\limits_{j=1}^{N_i}{n_{ij}}}{N_i+1})$ You need to normalize this formula to sum to 1 over the range of relevant Y. To get the country-wide estimate, you would need to use Monte-Carlo simulation for the sum of the $Y_i$ from each area based on the above formula. There are a couple inexpensive/free simulators out there.