#!/usr/bin/pythonimport sys,osimport cgiimport subprocessfrom subprocess import PIPE,Popenimport socketimport urllib2,urllibimport tempfileimport refrom PIL import Image, ImageDraw, ImageFontfaked_useragents ={ "Mozilla_unix_64": "Mozilla/5.0 (X11; U; Linux x86_64; de; rv:1.9.0.10) Gecko/2009042805 a Firefox 2.0.0.12 a Unix", "Mozilla_windows_xp": "Mozilla/5.0 (Windows; U; Windows NT 5.1; de; rv:1.9.0.10) Gecko/2009042316 Firefox/3.0.10 (.NET CLR 3.5.30729)", "IExplore_windows_xp": "Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 5.1; .NET CLR 2.0.50727; .NET CLR 3.0.4506.2152; .NET CLR 3.5.30729)"}# we need some default values# default_server:# key[0]: wms connect string# key[1]: projection system# key[2]: WMS/WEB --> are we querriing a webserver or an wms server# key[3]: dictonary with an cropbox# key[4]: dictonary with mapping from easy layernames -> list [real_layer_name,style_for layer]# key[5]: which user_agent should be faked# key[6]: which image format should be requesteddefault_server = {# whole germany - sometimes black"DOPV": ["http://gdz-dop.bkg.bund.de/arcgis/services/DOP/ImageServer/WMSServer?", "EPSG:25832", "WMS", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["0","default"], }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# nice orthos, but not every where these are colored "TH": ["http://www.geoproxy.geoportal-th.de/geoproxy/services?", "EPSG:31468", "WMS", {"top":0,"right":0,"bottom":20,"left":0}, { "DOPC":["dop20c","default"], "DOPSW":["dop20sw","default"], "DOP2C":["dop2c","default"], "DOP2SW":["dop2sw","default"], "GEMEINDEGRENZEN250":["GGK250","default"], "GEMEINDEGRENZEN":["gemeindegrenzen","default"], "GEMARKUNGSGRENZEN":["gemarkungsgrenten","default"], "FLURGRENZEN":["flurgrenzen","default"] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# sometimes the server don't delivers any answer, even if a tcp connection is# up and running"SN": ["http://atlas.egov.sachsen.de/geobak/servlet/gtEntryPoint;?", "EPSG:31468", "WMS", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["layer_3533",""], "UEBERSICHTSKARTE":["map_602",""], "STRASSEN":["layer_4419",""], "TOP":["layer_3866",""] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# when seeing these beach orthos ... ;-)"MV": ["http://www.gaia-mv.de/dienste/DOP?", "EPSG:31468", "WMS", {"top":5,"right":0,"bottom":0,"left":0}, { "DOPC":["DOP",""] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# yes we take the orthophotos from the web :-)"BY": ["http://www.geodaten.bayern.de/dyn/view/img", "EPSG:31468", "WEB", {"top":5,"right":0,"bottom":20,"left":0}, { "DOPC":["dop",""], "TK":["tk",""] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# oh dear, what a lousy webserver to query. the java version of theire client is also not working# until further explorations this server is commented out"SA": ["http://www.lvermgeo.sachsen-anhalt.de//de/geoservice/viewer/htmlviewer/htmlviewer.htm?powerProxyRequestURL=/navigator/basicviewer.jsp?", "EPSG:31468", "WEB", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["dop",""] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# when the server is up, it works well - but delivers only 40cm resolution"HE": ["http://thunder2.hessen.de/cgi-bin/lika-services/de-viewer/access/ogc-free-images.plx?", "EPSG:31467", "WMS", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["DOP04F1","default"] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# give me some working epsg-definition in /usr/share/proj/epsg and # the whole thing will be good"BB": ["http://isk.geobasis-bb.de/ows/dop20c.php?", "EPSG:325833", "WMS", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["18000",""] }, faked_useragents["Mozilla_unix_64"], "image/png"],# only here wile 'cause i checked the terrible 2m resolution :-D"BY2": ["http://www.geodaten.bayern.de/ogc/getogc.cgi?", "EPSG:4326", "WMS", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPC":["DOP",""] }, faked_useragents["Mozilla_unix_64"], "image/jpeg"],# not implemented yet"BW": ["http://rips-uis.lubw.baden-wuerttemberg.de/ims/cgi?ServiceName=badekarte&ClientVersion=4.0&Form=True&Encode=False", "EPSG:31467", "WEB", {"top":0,"right":0,"bottom":0,"left":0}, { "DOPBW":["dopbw",""] }, faked_useragents["Mozilla_unix_64"], "image/gif"]}default_wms_values = {"SERVICE": "WMS","VERSION": "1.1.1","REQUEST": "GetMap","BGCOLOR": "0xFFFFFF","TRANSPARENT": "TRUE","STYLES": "default","EXCEPTIONS": "application/vnd.ogc.se_inimage"}default_srs = "EPSG:4326"default_format ="image/png"default_precision = 8default_overlap = 4default_crop = Truedefault_fake_useragent = Truedefault_timeout = 30default_errorfont_name = "/usr/share/fonts/truetype/freefont/FreeMono.ttf"default_errorfont_size = 18# initialize variablesserver = None #Server to querry for mapwms = None #connect string for the remote wmsbbox = None #bbox paramet in llx lly urx ury as stringlayers = None #layers to request --> default_server[x][4]styles = None #styles for the layersprecision = None #precision in fractal when passing arguments srs_in = None #input spatial systemsrs_out = None #output spatial systemwidth = None #width height = None #heightoverlap = None #amount of enlargement the map in %dpsx = 0 #x-enlargement in srs_indpsy = 0 #y-enlargement in srs_ingrid_i = [] #input grid in srs_ingrid_o = [] #output grid in srs_ingrid_ot = [] #outpur grid in srs_outcrop = None #should the cropbox be recognizedcropbox = [] #cropbbox for cutting the imagefake_useragent = None #should the user-agent string be fakedtimeout = None #timeout in seconds we are waiting for an server responseformat_in = None #which image format should be deliveredformat_out = None #which format we are requesting# print the error strings#############def printerror (errorstring): if (width is None) or (height is None): print "Content-type: text/html" print "" print "<htlm><head><title>WMS-Errorfile</title></head><body>" print errorstring print "</body></html>" else: print "Content-type: image/jpeg\n" img = Image.new("RGB",(width,height)) if os.path.isfile(default_errorfont_name): font = ImageFont.truetype(default_errorfont_name,default_errorfont_size) else: font = ImageFont.load_default() draw = ImageDraw.Draw(img) for i,line in enumerate(errorstring.splitlines()): draw.text((5,15*i+5),line,(200,200,200),font) img.save(sys.stdout,"JPEG") sys.exit(2)# build an function for handeling an url request#############def fetchurl (req): if server in ["DOPV"]: #WMS Server from bund need those special headers req.add_header('Accept','text/*, image/png, image/gif, image/jpeg, *; q=.2, */*; q=.3') req.add_header('gdz-group','1028196') if fake_useragent: req.add_header('User-Agent',default_server[server][5]) try: answer = urllib2.urlopen(req) return answer except urllib2.URLError, e: printerror("Problem occured contacting outging web server %s\n%s" %(e,req.get_full_url())) except socket.error: printerror("The Server is not reachable")# transformes point ix,iy form srs_i to srs_o with precsision prec#############def transform (prec,srs_i, srs_o, ix,iy): param = '-s_srs %s -t_srs %s << START\n%.*f %.*f\nSTART' % (srs_i,srs_o,prec,ix,prec,iy) return subprocess.Popen(['gdaltransform '+param],stdout=PIPE,shell=True).communicate()[0].split(' ')[0:2]# search for necessary tools#############def tools_okay(): tools = ["gdaltransform","gdal_translate","gdalwarp"] missing = "" for tool in tools: try: subprocess.check_call(["which",tool]) except CalledProcessError: missing+="%s\n"% (tool) if missing != "": printerror("The following tools could not be found,\nor are not accessible from the cgi-script\n"+missing)############# parse the http GET parameters###########myparams = {}formdata = cgi.FieldStorage()for i in formdata: myparams[i.upper()]=formdata.getlist(i)for k in myparams: if k == "BBOX": bbox = myparams[k][0] if k == "SERVER": server = myparams[k][0] elif k == "LAYERS": layers = myparams[k][0] elif k == "STYLES": styles = myparams[k][0] elif k == "SRS": srs_in = myparams[k][0] elif k == "HEIGHT": height = myparams[k][0] elif k == "WIDTH": width = myparams[k][0] elif k == "CROP": crop = myparams[k][0] elif k in default_wms_values: default_wms_values[k] = myparams[k][0] elif k =="FAKEUSERAGENT": fake_useragent = myparams[k][0] elif k =="TIMEOUT": timeout = myparams[k][0] elif k =="FORMAT": format_in = myparams[k][0]# check the parameters##################if width is None: #check width printerror("There must be an width parameter")else: try: width = int(width) except ValueError: printerror("WIDTH must be an integer value")if height is None: #check height printerror("There must be an height parameter")else: try: height = int(height) except ValueError: printerror("HEIGHT must be an integer value")if bbox is None: #check bbox printerror("There must be an bbox field parameter") sys.exit(2)else: try: for i in bbox.split(','): grid_i.append(float(i)) except ValueError: printerror("There must be float values in bbox")if len(grid_i) != 4: #check iv ve have exact 4 values in bbox printerror("There must be 4 values in bbox")if server is None: #check for server from default_server printerror("There must be an SERVER parameter")else: server = server.upper() if server not in default_server: errorstr = "SERVER %s is unknow to the script.\n" % (server) errorstr += "Possible Values for Server are:\n" for s in default_server: errorstr += "%s\n" % (s) printerror(errorstr) else: #fill some variables with values form default_server wms = default_server[server][0] srs_out = default_server[server][1] format_out = default_server[server][6].split('/')[1]if layers is None: #check layer printerror("There must be an LAYERS parameter")else: reallayers = "" for l in layers.split(','): if l.upper() not in default_server[server][4]: errorstr = "LAYER '%s' is not defined\n" % (l) errorstr += "Possible Values for LAYERS are:\n" for l in default_server[server][4]: errorstr += "%s\n" %l printerror(errorstr) else: reallayers += default_server[server][4][l.upper()][0] givenlayers = layers #need this copy later for checking styles layers = reallayersif styles is None: #check styles realstyles = '' for l in givenlayers.upper().split(','): realstyles += default_server[server][4][l][1] styles = realstyleselse: if len(styles.split(',')) != len(layers.split(',')): printerror("The number of layers and styles in your request don't mach")if precision is None: #check precision precision = default_precisionif srs_in is None: #cehck srs srs_in = default_srsif width is None: #check width printerror("There must be an width parameter")else: try: width = int(width) except ValueError: printerror("WIDTH must be an integer value")if height is None: #check height printerror("There must be an height parameter")else: try: height = int(height) except ValueError: printerror("HEIGHT must be an integer value")if overlap is None: #check overlap overlap = default_overlapif crop is None: #check crop crop = default_cropelse: if crop.upper() == "TRUE": crop = True elif crop.upper() == "FALSE": crop = False else: crop = default_cropif fake_useragent is None: #check if we should fake the user agent fake_useragent = default_fake_useragentelse: if fake_useragent.upper() == "TRUE": fake_useragent = True elif fake_useragent.upper() == "FALSE": fake_useragent = False else: printerror("Value for FAKEUSERAGENT can only be True/False")if timeout is None: #check and set the timeout for socket operation timeout = default_timeoutelse: try: timeout = int(timeout) except ValueError: printerror("Vale for TIMEOUT can only be an integer")socket.setdefaulttimeout(timeout)if format_in is None: # check input graphics format format_in = default_format.split('/')[1]else: if re.match("image/.+",format_in) is not None: format_in = format_in.split('/')[1]print srs_inprint srs_out# now let's do the real work ;-)############################################################################################################ compute the enlargements for the 2 directionspsx = (grid_i[2]-grid_i[0])/width psy = (grid_i[3]-grid_i[1])/height# compute cropbox values for enlargement the aera and later cuttingif crop is True: values = default_server[server][3] cropbox = [psx*values["left"],psy*values["bottom"],psx*values["right"],psy*values["top"]]else: cropbox = [0,0,0,0] # compute the overlap values dpsx = (grid_i[2]-grid_i[0])/100*overlapdpsy = (grid_i[3]-grid_i[1])/100*overlap# compute the enlargment aera in original srs grid_o = [grid_i[0]-dpsx-cropbox[0],grid_i[1]-dpsy-cropbox[1],grid_i[2]+dpsx+cropbox[2],grid_i[3]+dpsy+cropbox[3]]# transform enlarged aerea into output srsgrid_ot = transform(precision,srs_in,srs_out,grid_o[0],grid_o[1])+transform(precision,srs_in,srs_out,grid_o[2],grid_o[3])#build the url for requesting the imageif default_server[server][2] == "WMS": add2url = "" for key in default_wms_values: add2url+="&%s=%s" % (key,default_wms_values[key]) add2url+="&FORMAT=image/%s&LAYERS=%s&BBOX=%s&SRS=%s&WIDTH=%i&HEIGHT=%i" % (format_out,layers,','.join(map(str,grid_ot)),srs_out,width,height) # get the image request = urllib2.Request(url=str(wms+add2url)) mapimage = fetchurl(request)# we need some special handling if we are getting the image from the webelse: if server in ["BY"]: request = urllib2.Request(url=str(wms)) #bavarian website uses only integer values, so we do round grid_otr = [float(x) for x in grid_ot] grid_otr = [int(grid_otr[0]),int(grid_otr[1]),int(grid_otr[2]+1),int(grid_otr[3]+1)] #add the request parameters 2 the post data request.add_data("url=%s_%s_%s_%s_%s" % (layers,','.join(map(str,grid_ot)),width,height,format_out)) picturepath = fetchurl(request) #get the picture request2 = urllib2.Request(url="%s://%s%s" % (request.get_type(),request.get_host(),picturepath.read())) mapimage = fetchurl(request2) else: errorwrite("Request of these webserver\n is still not implemented! Sorry")# go for image processing (get, crop, warp, translate etc.########################################################################################################wmsinput = tempfile.NamedTemporaryFile(mode='w+b',suffix='.%s.tmp'%(format_out),prefix='wmsviewer-')wmsinput.write(mapimage.read())wmsinput.flush()wmsinput.seek(0)if wmsinput.read(1) == "": printerror("The outgoing server did not provide any data within time")#create some temporary file for image translationgdaltrans1 = tempfile.NamedTemporaryFile(mode='w+b',suffix='.tif.tmp',prefix='wmsviewer-')gdalwarp = tempfile.NamedTemporaryFile(mode='w+b',suffix='.tif.tmp',prefix='wmsviewer-')gdaltrans2 = tempfile.NamedTemporaryFile(mode='w+b',suffix='.%s.tmp'%(format_in),prefix='wmsviewer-')#build a georeferenced imageparamgdaltrans1 = " -a_srs %s -a_ullr %s %s %s %s %s %s" % (srs_out,grid_ot[0],grid_ot[3],grid_ot[2],grid_ot[1],wmsinput.name,gdaltrans1.name)try: returncode = subprocess.call('gdal_translate'+paramgdaltrans1,shell=True,stdout=PIPE) if returncode != 0: printerror("Building the georeferenced image failed with returncode :%i"% (returncode))except OSError,e: errorstr = "Building the georeferenced image gailed\n" errorstr += "Error:%s\n" errorstr += "process: gdal_translate\n" errorstr += "args:%s\n" % (e,paramgdaltrans1) printerror (errorstr)gdaltrans1.flush()#translate to srs_inparamgdalwarp = " -q -s_srs %s -t_srs %s %s %s" % (srs_out,srs_in,gdaltrans1.name,gdalwarp.name)subprocess.call('gdalwarp'+paramgdalwarp,shell=True,stdout=PIPE)gdalwarp.flush()#crop and translate to output formatparamgdaltrans2 =" -quiet -outsize %i %i -of %s -projwin %f %f %f %f %s %s" % (width,height,format_in,grid_i[0],grid_i[3],grid_i[2],grid_i[1],gdalwarp.name,gdaltrans2.name)proc = subprocess.call('gdal_translate'+paramgdaltrans2,shell=True,stdout=PIPE)gdaltrans2.flush()gdaltrans2.seek(0)#deliver the image to the clientprint "Content-Type: image/%s\n" % (format_in)print gdaltrans2.read()#do some cleanup ;-)if os.path.isfile("%s.aux.xml"%(gdaltrans2.name)): os.unlink("%s.aux.xml"%(gdaltrans2.name))#python could do this himself, but better safe, than sorrygdaltrans2.close()gdalwarp.close()gdaltrans1.close()wmsinput.close()