function get_hf_data, which_layer, bounding_box, outfile=fname, wcs=WCS, wms=WMS, no_display=NO_DISPLAY, get_layers=GET_LAYERS ; this function reads a wms_layer for a given url for the ; get_hf_data( [layer number], [bounding box], [outfile=fname], /WCS, /WMS, /GET_LAYERS) ; ; How to use this program: ; ; to determine what layers there are, you can type the following at ; the IDL prompt: ; ; IDL> print, get_hf_data() ; IDL> layer_wcs = get_hf_data(/wcs,/get_layers) ; ; to find the bounding box data for any of the layers, you type: ; ; IDL> bb = get_hf_data(0, /wcs) ; for the WCS layer number zero ; ; you can edit the bounding box layer to narrow down the region you ; are looking at ; ; IDL> print, bb ; IDL> bb = bb[1] ; this is a UTM layer, but you can just as easily edit the lat/lon layer ; IDL> bb.minx = 730000 ; IDL> bb.maxx = 732000 ; IDL> bb.miny = 4712000 ; IDL> bb.maxy = 4714000 ; ; or IDL> bb.minx = 730000 & bb.maxx = 732000 & bb.miny = 4712000 & bb.maxy = 4714000 ; ; and now, you can create the image file ; ; IDL> bbimg = get_hf_data(0, bb, outfile = 'canopy_height.tiff', /wcs) ; IDL> tvscl, bbimg.image ; ; This program is written by Paul Siqueira, University of Massachusetts, Amherst return_value = ' ' if ( (not Keyword_set(WCS)) AND (not Keyword_set(WMS)) )then begin WCS = 1 ; the default is to return information from both WCS and WMS WMS = 1 endif fname_default = 'temp_file' if(n_elements(fname) eq 0) then fname = fname_default ; if the Harvard Forest URL changes, then this next line should be changed url_wcs = 'http://107.20.174.13/ArcGIS/services/Harvard_Forest_Remote_Sensing/MapServer/WCSServer?' url_wms = 'http://107.20.174.13/ArcGIS/services/Harvard_Forest_Remote_Sensing/MapServer/WMSServer?' owcs = obj_new('idlnetogcwcs',COVERAGE_FILENAME=fname) ; create an object for reading maps owms = obj_new('idlnetogcwms') owms->ParseURL, url_wms owcs->ParseURL, url_wcs number_layers_wcs = owcs->GetCapabilities() number_layers_wms = owms->GetCapabilities() layer_wms = owms->GetLayers() layer_wcs = owcs->GetCoverageOfferingBriefs() spaces = ' ' ; if no layer number is given, provide a review of all of the ; available layers if (n_elements(which_layer) eq 0) then begin if(Keyword_Set(WMS))then begin print, ' ' print, 'WMS Layers' print, '--------------------------------' for ii = 0, number_layers_wms-1 do begin print, strtrim(ii,2) + ' ' + strmid([layer_wms[ii].title+spaces],0,20) + ' ' + 'xxx' endfor print, ' ' endif if(Keyword_Set(WCS))then begin print, 'WCS Layers' print, '--------------------------------' for ii = 0, number_layers_wcs-1 do begin print, strtrim(ii,2) + ' ' + strmid([layer_wcs[ii].label+spaces],0,20) + ' ' + 'xxx' endfor endif if (keyword_set(GET_LAYERS)) then begin if(keyword_set(WCS)) then return_value = {WCS:layer_wcs} if(keyword_set(WMS)) then return_value = {WCS:layer_wms} if(keyword_set(WCS) AND keyword_set(WMS)) then return_value = {WCS:layer_wcs, WMS:layer_wms} endif else begin return_value = ' ' endelse endif else begin ; a layer number (which_layer) has been provided ll_fmt = '(%"%s: %8.4f %8.4f")' meters_fmt = '(%"%s: %d %d")' if(keyword_set(WMS))then begin nbbox = layer_wms[which_layer].num_bounding_box bbox = layer_wms[which_layer].Bounding_Box tags = tag_names(bbox) print, ' ' print, ['Layer name: ' + layer_wms[which_layer].title] for ibox = 0, nbbox-1 do begin if (ibox eq 0) then begin print, format = '(%"\n%s")', bbox[ibox].srs print, format = ll_fmt, 'min: ',bbox[ibox].minx, bbox[ibox].miny print, format = ll_fmt, 'max: ',bbox[ibox].maxx, bbox[ibox].maxy endif else begin print, format = '(%"\n%s")', bbox[ibox].srs print, format = meters_fmt, 'min: ',round(float(bbox[ibox].minx)), round(float(bbox[ibox].miny)) print, format = meters_fmt, 'max: ',round(float(bbox[ibox].maxx)), round(float(bbox[ibox].maxy)) endelse print, ' ' endfor endif else begin ; for the WCS keyword instead (normally executed) min_vals = strsplit(layer_wcs[which_layer].pos1,/EXTRACT) max_vals = strsplit(layer_wcs[which_layer].pos2,/EXTRACT) minx = float(min_vals[0]) miny = float(min_vals[1]) maxx = float(max_vals[0]) maxy = float(max_vals[1]) utm = ll_to_utm( [minx,maxx],[miny,maxy], datum=2) utm_str = string(utm) utmzone = cgutmzone((minx + maxx)/2.0,(miny+maxy)/2.0,formal=str_designation) bbox = {SRS:'EPSG:4326', CRS:' ', MINX:min_vals[0], MINY:min_vals[1], MAXX:max_vals[0], MAXY:max_vals[1], RESX:' ', RESY:' '} bbox = [bbox, bbox] bbox[1] = {SRS:'UTM:18T', CRS:' ', MINX:utm_str[0,0], MINY:utm_str[0,1], MAXX:utm_str[1,0], MAXY:utm_str[1,1], RESX:' ', RESY:' '} print, ' ' print, 'Latitude/Longitude coordinates' print, '------------------------------' print, format = ll_fmt, 'min: ',minx, miny print, format = ll_fmt, 'max: ',maxx, maxy print, ' ' print, ['UTM coordinates (zone=' + str_designation + ')'] print, '------------------------------' print, format = meters_fmt, 'min: ',utm[0,0],utm[1,0] print, format = meters_fmt, 'max: ',utm[0,1],utm[1,1] endelse return_value = bbox srs_wgs84 = 'EPSG:4326' if (n_elements(bounding_box) ne 0) then begin srs = bounding_box[0].srs minx = double(bounding_box[0].minx) miny = double(bounding_box[0].miny) maxx = double(bounding_box[0].maxx) maxy = double(bounding_box[0].maxy) if(strmatch('UTM:18T',strupcase(srs))) then begin east = [minx, maxx] north = [miny,maxy] ll = utm_to_ll(east,north,'WGS84',zone=18) minx = ll[0,0] miny = ll[1,0] maxx = ll[0,1] maxy = ll[1,1] endif if(minx gt maxx) then begin temp = minx minx = maxx maxx = temp endif if(miny gt maxy) then begin temp = miny miny = maxy maxy = temp endif utm = ll_to_utm( [minx,maxx],[miny,maxy], datum=2) minx = string(minx) maxx = string(maxx) miny = string(miny) maxy = string(maxy) bb = {srs:srs, minx:minx, miny:miny, maxx:maxx, maxy:maxy} print,' ' print,'New values of the bounding box' print,'==============================' print,'Latitude/Longitude coordinates' print,'------------------------------' print, format = ll_fmt, 'min: ',minx, miny print, format = ll_fmt, 'max: ',maxx, maxy print, ' ' print, ['UTM coordinates (zone=' + str_designation + ')'] print, '------------------------------' print, format = meters_fmt, 'min: ',utm[0,0],utm[1,0] print, format = meters_fmt, 'max: ',utm[0,1],utm[1,1] query_bbox = bb.minx + ',' + bb.miny + ',' query_bbox = query_bbox + bb.maxx + ',' + bb.maxy query_bbox = strcompress(query_bbox,/REMOVE_ALL) if(n_elements(width) eq 0) then width = 512 if(n_elements(height) eq 0) then height = 512 if(keyword_set(WMS))then begin query_string = 'LAYERS=' + layer_wms[which_layer].name query_string = query_string + '&STYLES=' + layer_wms[which_layer].style.name query_string = query_string + '&SRS=' + srs_wgs84 query_string = query_string + '&BBOX=' + query_bbox query_string = query_string + '&WIDTH=' + string(width) query_string = query_string + '&HEIGHT=' + string(height) query_string = query_string + '&FORMAT=' + layer_wms[which_layer].map_format[2] ; this is 'image/tiff' map_file = owms->GetMap(query_string) print, 'file retrieved = ', map_file endif else begin image_format = 'GeoTIFF' image_resx = 1.0 image_resy = 1.0 crs_identifier = srs_wgs84 query_string = 'Coverage=' + layer_wcs[which_layer].name query_string = query_string + '&CRS=' + crs_identifier query_string = query_string + '&BBOX=' + query_bbox query_string = query_string + '&FORMAT=' + image_format query_string = query_string + '&WIDTH=' + strtrim(string(width),2) query_string = query_string + '&HEIGHT=' + strtrim(string(height),2) map_file = owcs->GetCoverage(query_string) print, 'file retrieved = ', map_file img = read_tiff(map_file,geotiff=geotiff_struct) return_value = {image:reverse(img,2), geotiff:geotiff_struct} if(not Keyword_set(NO_DISPLAY)) then begin wtitle = ['Window ' + strcompress(string(!D.window),/remove_all) + ':' + layer_wcs[which_layer].label] window, !D.window, xsize=width, ysize=height, title = wtitle tvscl, img endif endelse endif endelse if(n_elements(return_value) eq 0) then return_value = number_layers return, return_value end