Calculate the scale given a bounding box and the width and the height of a image

I need to calculate the scale given a bounding box in map coordinates (EPSG:900913) and the width/height of the image.

This is becose Im implementing a very basic WMS server but in certain moments i need to know the scale of the image.

A simple request would be :

  wms.exe?LAYERS=all&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fjpeg&SRS=EPSG%3A900913&BBOX=-7810029.800979,-2298920.0623849,-7809724.0528659,-2298614.3142718&WIDTH=256&HEIGHT=256

that request is originated from a openlayer component in a web page configured with this parameters :

var optionsGmap = {
projection: new OpenLayers.Projection("EPSG:900913"),
units: "m",
sphericalMercator: true,
numZoomLevels: 16,
maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
scales : [110936068.185464,55468034.0927321,27734017.0463661,13867008.523183,6933504.26159152,
    3466752.13079576,1733376.06539788,866688.03269894,433344.01634947,216672.008174735,108336.004087367,54168.0020436837,27084.0010218419,
    13542.0005109209,6771.00025546046,3385.50012773023],
maxResolution: "auto" };

theres formula to calculate this? I have access to a postgres+posgis in the backend but I prefer to do the calculation by my self in C code.

Solution

The lengths of each side of an orthogonal geographical area are as follows:

lengthEastWestNorth = (east - west) * ((equatorial_circumference * cos(north)) / (2*pi) ) ;
lengthEastWestSouth = (east - west) * ((equatorial_circumference* cos(south)) / (2*pi)) ;
lengthNorthSouth = (north - south) * (polar_circumference/(2*pi)) ;

Note that east, west, north, south are long/lat angles in radians since that is what the trig function in C's <math.h> require. Conversion from degrees to radians is simply r = pi/180.0. The units of the result are whatever units you use for the circumferences. The north-south sides are of equal length, the east-west dimensions differ depending on latitude, and the length of the north-south extent, hence three values not two or four.

To calculate the scale you probably need only one of these values, I suggest the lengthNorthSouth value since it requires no trig functions and therefore no conversion to radians (divide by 360 rather than 2*pi in this case).

The length of one side divided by the corresponding side length of the image will yield the scale.

Geographical calculations are necessarily approximations, and the map projection may use different circumferential values for example, so you will have to determine the nearest preferred scale to determine the intended scale. Subtract the calculated scale from each of the preferred scales in turn to determine the scale with the smallest absolute (remove the sign) difference.