function get_massgis_data, which_layer, bbox=bbox_in, resolution=resolution, all_layers=ALL_LAYERS, get_data=GET_DATA, filename=fname, no_display=NO_DISPLAY ; function get_massgis_data ; ; usage: print, get_massgis_data() ; print, get_massgis_data(/all_layers) ; layer = get_massgis_data(0) ; ; bb_shaler = create_bounding_box(-72.192,42.530,2000,1500) ; bb_prosp = create_bounding_box(-72.193,42.519,3000,4000) ; bb_simes = create_bounding_box(-72.220,42.463,2000,4000) ; data = get_massgis_data(0,bbox=bb) ; data = get_massgis_data(0,bbox=bb,resolution=1) ; data = get_massgis_data(0,resolution=0.5) ; data = get_massgis_data(0,/get_data) ; data = get_massgis_data(0,/get_data,filename='ortho_2001') ; data = get_massgis_data(0,/get_data,bbox=bb,resolution=1,filename='ortho_2001') ; ; The input variables have the following meanings ; ; get_massgis_data(which_layer,bbox=BB,resolution=res,filename=fname,/all_layers,/get_data,/no_display) ; ; which_layer - refers to the numbers that are in the ; first column that shows up when the user ; types: print, get_massgis_data() ; bbox - a bounding box ; resolution - the resolution of the returned data (units in meters) ; /all_layers - a keyword, when set, returns all of the ; layers. To be used with the print, ; get_massgis_data(/all_layers) command. This ; keyword does not retrieve actual data. ; /get_data - normally, ; /no_display - disables the display of the retrieved data ; ; ; ; ; ; ; ; This program is written by Paul Siqueira, University of Massachusetts, Amherst if (n_elements(bbox_in) ne 0) then begin bbox = bbox_in endif fname_default = 'MassGIS_file' if(n_elements(fname) eq 0) then fname = fname_default owcs = obj_new('idlnetogcwcs',COVERAGE_FILENAME=fname) ; create an object for reading maps url = 'http://giswebservices.massgis.state.ma.us/geoserver/wcs/WCSServer?' owcs->ParseURL, url nlayers = owcs->GetCapabilities() layer_wcs = owcs->GetCoverageOfferingBriefs() ; print, tag_names(layer_wcs) if (n_elements(which_layer) eq 0) then begin if(keyword_set(ALL_LAYERS))then begin for ilayer = 0, nlayers-1 do begin outstring = strtrim(ilayer,2) + ' ' + layer_wcs[ilayer].label print, outstring endfor endif else begin interesting_layers = [0,1,3,5,12,95,96,97,102,103,104,105,106] ; these are the interesting layer numbers and names as of July 2013 ; 0 2001 Color Orthos ; 1 2005 Color Orthos ; 3 2009 Color Orthos 30cm ; 5 Black and White Orthos (1990s) ; 12 Elevation ; 95 Shaded Relief ; 96 Shaded Relief 2005 ; 97 Shaded Relief New England ; 102 USGS Topos ; 103 Wind Speed at 100M ; 104 Wind Speed at 30M ; 105 Wind Speed at 50M ; 106 Wind Speed at 70M for ii = 0, n_elements(interesting_layers)-1 do begin ilayer = interesting_layers[ii] outstring = layer_wcs[ilayer].index + ' ' + layer_wcs[ilayer].label print, outstring endfor endelse return_value = ' ' endif else begin ; in other words, there is a layer number specified (which_layer) return_value = layer_wcs[which_layer] get_data_flag = 0 if(n_elements(bbox) ne 0) then get_data_flag = 1 if(n_elements(resolution) ne 0) then get_data_flag = 1 if(keyword_set(GET_DATA)) then get_data_flag = 1 if(get_data_flag eq 1) then begin if(n_elements(bbox_in) eq 0) then begin default_minx = -72.193390 default_miny = 42.528343 default_minx = -72.19100 default_miny = 42.53100 if(n_elements(resolution) eq 0) then begin default_deltax = 1000 default_deltay = 1000 bbox = create_bounding_box(default_minx, default_miny, default_deltax, default_deltay) bbox = bbox[0] endif else begin ; resolution is set, but bbox is not set default_deltax = 512*resolution default_deltay = 512*resolution bbox = create_bounding_box(default_minx, default_miny, default_deltax, default_deltay) bbox = bbox[0] endelse endif if(n_elements(resolution) eq 0) then begin width = 512 height = 512 endif else begin ; the resolution is given in meters, so we need to convert that into lat/lon utm = ll_to_utm([bbox[0].minx, bbox[0].maxx],[bbox[0].miny, bbox[0].maxy],datum=2) deltax = utm[0,1] - utm[0,0] deltay = utm[1,1] - utm[1,0] ; print, 'deltax: ', deltax ; print, 'deltay: ', deltay width = round(deltax/resolution) height = round(deltay/resolution) endelse bbox = bbox[0] minx = string(bbox.minx) miny = string(bbox.miny) maxx = string(bbox.maxx) maxy = string(bbox.maxy) query_bbox = minx + ',' + miny + ',' query_bbox = query_bbox + maxx + ',' + maxy query_bbox = strcompress(query_bbox,/REMOVE_ALL) image_format = 'GeoTIFF' crs_identifier = layer_wcs[which_layer].srs_name crs_identifier = 'EPSG:4326' crs_identifier = 'CRS:84' 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) img = read_tiff(map_file,geotiff=geotiff_struct) ; return_value = {image:reverse(img,2), geotiff:geotiff_struct} if (strcmp(fname,fname_default)) then begin file_delete, map_file endif else begin print, 'file retrieved = ', map_file endelse dims = size(img) if(dims[0] eq 3)then begin img = reverse(img,3) endif else begin img = reverse(img,2) endelse return_value = {image:img, geotiff:geotiff_struct} if(not Keyword_set(NO_DISPLAY)) then begin if(!D.window eq -1) then begin wtitle = ['Window ' + strcompress(string(0),/remove_all) + ':' + layer_wcs[which_layer].label] window, xsize=width, ysize=height, title = wtitle endif else begin wtitle = ['Window ' + strcompress(string(!D.window),/remove_all) + ':' + layer_wcs[which_layer].label] window, !D.window, xsize=width, ysize=height, title = wtitle endelse if(dims[0] eq 3)then begin tvscl, return_value.image, true=1 endif else begin tvscl, return_value.image, /nan endelse endif endif endelse return, return_value end