Calculating the volume of a restaurant take-away box that is circular on the bottom and square on the top
A different approximation could be to take the cross-section at each height not as a linear interpolation between top and bottom surface, but as squares with rounded corners. This fits the photograph more closely, as the linear interpolation approach has a discontinuity in the radius of curvature at the bottom: The point that at the top is an edge of the square is modelled a crease along the entire side, while in the photograph, no such crease can be seen.
If we take the cross-section to have the above shape of a rounded square, we can calculate its area by using the formulas for circles and rectangles:
$$A = \pi \cdot r^2 + 2 \cdot a \cdot b - b^2$$
Using $b=a-2\cdot r$, we can eliminate that variable:
$$A = a^2 + (\pi -4) \cdot r^2$$
Now we assume that those variables vary linearly between top and bottom: $r$ goes from $r_0$ at the bottom to $0$ at the top, and $a$ goes from $2r_0$ at the bottom to $a_0$ at the top, giving us:
$$r(h) = \left(1-\frac{h}{h_0}\right)\cdot r_0$$
$$a(h) = \left(1-\frac{h}{h_0}\right)\cdot 2r_0 + \frac{h}{h_0} \cdot a_0$$
and thus:
$$A(h) = \left(1-\frac{h}{h_0}\right)^2 \cdot \pi \cdot r_0^2 + \left(1-\frac{h}{h_0}\right) \cdot \frac{h}{h_0} \cdot 4 \cdot r_0 \cdot a_0 + \left(\frac{h}{h_0}\right)^2 \cdot a_0^2$$
or, better ordered for integration:
$$A(h) = \pi \cdot r_0^2 + \frac{h}{h_0} \cdot \left( 4\cdot r_o \cdot a_0 - 2 \cdot \pi \cdot r_0^2 \right) + \left(\frac{h}{h_0}\right)^2 \cdot \left( a_0^2 - 4\cdot r_0 \cdot a_0 + \pi \cdot r_0^2\right) $$
Which gives us then:
$$ V = \int_0^{h_0} A(h) \mathrm{d}h= h_0 \cdot \pi \cdot r_0^2 + \frac{1}{2} \cdot \frac{h_0^2}{h_0} \cdot \left( 4\cdot r_o \cdot a_0 - 2 \cdot \pi \cdot r_0^2 \right) + \frac{1}{3}\cdot\frac{h_0^3}{h_0^2} \cdot \left( a_0^2 - 4\cdot r_0 \cdot a_0 + \pi \cdot r_0^2\right)$$
$$=h_0\cdot\left(\frac{1}{3}\pi r_0^2 + \frac{2}{3} r_0 a_0 + \frac{1}{3} a_0^2\right)$$
And here is a quick render of what this looks like in 3d:
This result can be obtained classically without any integrals as well, using only the formulas for the volume of conic solid and for cuboids. To do this, we split the solid into nine parts:
The central, yellow part is simply a square pyramid, and has a volume of $\frac{1}{3}\cdot h_0\cdot a_0^2$. The four magenta pieces are oblique cones with a quarter-circle as base. Again, using the formula for conic solids, their volume is each $\frac{1}{12}\cdot h_0\cdot\pi\cdot r_0^2$
And the four cyan pieces are irregularly shaped tetrahedrons, but we can determine their volume by adding a few pieces to see how they fit into a cuboid of the measurements $a_0\times r_0 \times h_0$:
In the more general case of $\frac{a_0}{2} \neq r_0$, there will be a fourth green piece needed (which would go in front and block our view of everything). However, the green pieces together form a prisma of height $a_0$ and a top surface area of $\frac{1}{2} \cdot r_0 \cdot h_0$, and the two blue-grey pieces are two pyramids with height $h_0$ and rectangular bases of size $\frac{1}{2} a_0 \times r_0$. Therefore, the volume of the tetrahedron is:
$$a_0 \cdot h_0 \cdot r_0 - a_0 \cdot \frac{1}{2} \cdot r_0 \cdot h_0 - 2 \cdot \frac{1}{3} \cdot h_0 \cdot \frac{1}{2} a_0 \cdot r_0 = \frac{1}{6} h_0 \cdot r_0 \cdot a_0$$
Putting it all together, we get the volume as:
$$V=\frac{1}{3}\cdot h_0\cdot a_0^2 + 4\cdot \frac{1}{12}\cdot h_0\cdot\pi\cdot r_0^2 + 4 \cdot \frac{1}{6} h_0 \cdot r_0 \cdot a_0$$ $$=h_0\cdot\left(\frac{1}{3}\pi r_0^2 + \frac{2}{3} r_0 a_0 + \frac{1}{3} a_0^2\right)$$
A more smooth model along the same lines would be to interpolate superellipses between the bottom and top, which similarly would give a creaseless change in the radius of curvature beneath the corners. However, superellipse areas have a formula involving the gamma function, and are therefore not easy to integrate again to get a volume.
As the shape of the solid is not clearly defined, I'll make the simplest assumption: lateral surface is made of lines, connecting every point $P$ of square base to the point $P'$ of circular base with the same longitude $\theta$ (see figure below).
In that case, if $r$ is the radius of the circular base, $2a$ the side of the square base, $h$ the distance between bases, a section of the solid at height $x$ (red line in the figure) is formed by points $M$ having a radial distance $OM=\rho(\theta)$ from the axis of the solid given by: $$ \rho(\theta)={a\over\cos\theta}{x\over h}+r\left(1-{x\over h}\right), \quad\text{for}\quad -{\pi\over4}\le\theta\le{\pi\over4}. $$ A quarter of the solid has then a volume given by: $$ {V\over4}=\int_0^h\int_{-\pi/4}^{\pi/4}\int_0^{\rho(\theta)}r\,dr\,d\theta\,dx= \frac{h}{12} \left(4 a^2+\pi r^2 + 2ar \ln{\sqrt2+1\over\sqrt2-1}\right). $$
Since it differs in circumference
Are you sure about this? One thing to note about the box is are the vertical stripes. Your image suggests that they have equal width and I see production- or design-wise reason why they shouldn’t have. However, this would make stacking them impossible (as noted by AlienAtSystem), if the boxes come partially pre-folded – which is hard to say. For this answer, I assume that cross sections of the box must have the same circumference, in particular the bottom and the top.
If the bottom is a circle with radius $r$, its area is $A_\text{bot} := πr^2$ and its circumference is $c := 2πr$. If the top is a square, it must have edges of length $a := \frac{c}{4} = \frac{πr}{2}$ and thus an area of $A_\text{top} := a^2 = \frac{π^2r^2}{4}.$
Now, for all the other horizontal cross sections, we have to make assumptions regarding their shape. The easiest assumption is that their area changes linearly along the height of the box¹. With this, we get that the average area of a horizontal cross section is $A_\text{av} = \frac{1}{2} (A_\text{top}+A_\text{bot})$. With $h$ as the height of the box, we arrive at:
$$V = A_\text{av}h = \frac{h}{2} (A_\text{top}+A_\text{bot}) = \frac{πhr^2 (π+4)}{8}.$$
Since there is little difference between the top and bottom area ($\frac{A_\text{bot}}{A_\text{top}} = \frac{4}{π} ≈ 1.27$), this approximation should suffice for most practical purposes.
If my initial assumption should be wrong, you can do the same calculation with $a$ decoupled from $c$ and arrive at $V=\frac{1}{2} h (πr^2+a^2)$.
¹ One might argue that it is the square root of the area should change linearly, as it does for cones, pyramids, and similar, but that is as much an assumption as a linearly changing area, given that we have bulging and similar.