1 / 46

Basic converting and combining data sets in v10.x

Basic converting and combining data sets in v10.x. USGS MRCTR Lab (Trent Hare), Jan. 2013 – GIS Workshop. Mars - first glance the data sets are overwhelming . Mariner, Viking (lander, orbiter), Hubble Mars Pathfinder: IMP, APXS, ASI/MET Mars Global Surveyor: MOC, TES, MOLA, MAG/ER

goldy
Download Presentation

Basic converting and combining data sets in v10.x

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Basic converting and combining data sets in v10.x USGS MRCTR Lab (Trent Hare), Jan. 2013 – GIS Workshop

  2. Mars - first glance the data sets are overwhelming... • Mariner, Viking (lander, orbiter), Hubble • Mars Pathfinder: IMP, APXS, ASI/MET • Mars Global Surveyor: MOC, TES, MOLA, MAG/ER • Mars Odyssey: THEMIS, GRS, NS, MARIE, HEND • Mars Express: HRSC, OMEGA, MARSIS/ASPERA/SPICAM • Mars Exploration Rovers: Pancam, Mini-TES, MB, APXS, MI, RAT, Navcam, Hazcam • Mars Reconnaissance Orbiter: HiRISE, CRISM, CTX, MARCI, MCS, SHARAD • Phoenix: RAD, MARDI, SSI, TEGA, MECA, MET • Mars Science Laboratory: MastCam, MARDI, MAHLI, ChemCam, APXS, CheMin, SAM, DAN, RAD, REMS

  3. How to use low-level PDS • Low-level PDS images (e.g. EDRs) are basically “raw” – no map projection – you should notbring it into a GIS • How do you map project EDR PDS images ISIS - Integrated Software for Imagers and Spectrometers Suse Linux, Solaris UNIX, Mac OSX -- (Linux Virtual Machine) http://isis.astrogeology.usgs.gov/ (2013 – new ISIS3 map projection web service) VICAR - Video Image Communication And Retrieval http://www-mipl.jpl.nasa.gov/ HRSC version maintained at DLR 3 GIS for Planetary Mappers

  4. RULES of the GIS ROAD • For ISIS processing • Best to set sameprojection and parameters for all • Note: optional to set same resolution • For visual (thematic) images, best to convert to 8bit • For “data” (e.g. DEM, Temperature -- 16,32 bit), set all ISIS Special Pixel Values to NULL (using specpix, stretch, bit2bit) • For global • If lonsys=360, then set clon=180 • If lonsys=180, then set clon=0 (better supported) • Don’t use funky projections

  5. Displaying 16, 32 bit ISIS cubes 1.) Right click layer, Select “properties” 2.) Symbology tab, select stretch: Std Dev. 3.) Yes to “calculate stats”, hit okay

  6. Batch Statistics open or drag files into here Tip: to calculate faster for large images change to 2, 10, 100, …

  7. When the range still looks bad - after stats

  8. First: try to set the NoDATA value in ArcCatalog Right click image, “Properties…” Type in value or try “Compute”

  9. Second: calculate valid range (setnull)

  10. Results for: calculate valid range (setnull)

  11. Or just Batch SET NULL for all images • USGS Image Toolbox ( http://bit.ly/q33Vqa ) • written for 9.x but works for 10

  12. No need for 16, 32? – convert to 8bit • ISIS • Isis2std - easiest but • Cannot convert files over 4GB (except Jpeg2000) • Tiff currently uses deflate compression – might need work • Can only convert to 8bit (except Jpeg2000 -8 or 16bit) • Does not support embedded projection (just worldfile) • Bit2bit (new @ ISIS 3.3.0) • Reason: convert from 32bit to 16 or 8bit ISIS cube • Requires app to still support ISIS format • Linear stretches only? • Mapping of clipped percentages to NULLs – generally not good

  13. No need for 16, 32? – convert to 8bit • ISIS – “Stretch” seems to be the best. After run, bring use ISIS cub in ArcMap. Stretch method 1>stretch from=input.cub to=output_8bit.cub+8bit+1:254 USEPERCENTAGES=true pairs="0:1 100:254" null=0 lis=0 lrs=0 his=255 hrs=255This allows you to specify input percentages for the mapping pairs. Thus when USEPERCENTAGES=true is set pairs="0:1 100:254" means:map 0% to 1 (or the file's min value to 1) and 100% to 254 (file's max value).Stretch method 2This also means you can apply a recommended 0.5% clip to remove the potential extraneous lows and highs like: > stretch from=input.cub to=output_8bit.cub+8bit+1:254 USEPERCENTAGES=true pairs="0:1 0.5:1 99.5:254 100:254" null=0 lis=0 lrs=0 his=255 hrs=255

  14. No need for 16, 32 – convert to 8bit • GDAL method • Manual method (more: http://bit.ly/pprlMK) > gdalinfo –mm in.cub (returns min/max, now convert) > gdal_translate –ot byte –scale min max 1 255 –a_nodata 0 in.cubout.tif > gdal_translate–of PNG –ot byte –scale min max 1 255 –a_nodata 0 in.cub out.png > gdal_translate–of JP2KAK –co quality=100 –ot byte –scale min max 1 255 –a_nodata 0 in.cubout.jp2 • On Astro Machines (more: http://bit.ly/oxIsQ7) > to8bit_gdal_tif.csh in.cubout.tif > to8bit_gdal_jp2.csh in.cub out.jp2 > to8bit_gdal_png.csh in.cubout.png

  15. GDAL for 32bit Map Projected PDS and ISIS3 GDAL (binaries available using FWtools and OSGeo4W): > gdal_translate –of GTIFF isis_ver3.cub isis_ver3.tif GDAL Tips: https://isis.astrogeology.usgs.gov/IsisSupport/index.php/topic,2172.0.html 15 GIS for Planetary Mappers

  16. Batch Command Line Tip Batch conversion Tips: Unix/Linux example: -------------------------------------------------------------------------------- foreachi (*.cub) foreach> perl isis2world.pl -e $i foreach> end In MsDOS command window loop (for Windows machines) example: -------------------------------------------------------------------------------- for %i in (*.cub) do isis2world -e %i 16 GIS for Planetary Mappers

  17. Trouble with ISIS cubes in ArcMap/GDAL? • Back-up (old) conversion method • First run isis2raw or isis2std (on “in.cub”) Now run • For raw run: > isis3world.pl –e –prjin.cub • For png run: > isis3world.pl –p –prjin.cub • For tiff run: > isis3world.pl –t –prjin.cub • For jpeg: > isis3world.pl –j –prjin.cub • For jpeg2000: > isis3world.pl –J –prjin.cub You will then need to assign the created projection to the output file using the new *.prj file. There are batch methods available: USGS Image Toolbox ( http://bit.ly/q33Vqa ) Link for script: Isis3world.pl

  18. Hands-on (HiRISE – Exercise 01) • Note: most HiRISE Jp2 images have embedded geospatial labels. However, they contain a Equirectangular parameter flaw when using GDAL. To fix, remap this parameter for ArcMap Open command prompt in directory with image and run: > fix_jp2 PSP_004365_1745_COLOR.JP2 Related: https://isis.astrogeology.usgs.gov/IsisSupport/index.php/topic,2339.0.html https://isis.astrogeology.usgs.gov/IsisSupport/index.php/topic,3440.msg13384.html • By setting the data frame to the projection of the HiRISE image you can now further nudge the HiRISE registration using ArcMap’s Georegistration Toolbar. 18 GIS for Planetary Mappers

  19. Simple Image RegistrationUsing a GIS Worldfile 19 GIS for Planetary Mappers

  20. * Worldfile • Most simple image registration 5.0 (size of pixel in x direction) – A 0.0 (rotation term for row) - D 0.0 (rotation term for column) - B -5.0 (size of pixel in y direction) - E 492169.690 (x coordinate of center of upper left pixel in map units) - C 54523.3180 (y coordinate of center of upper left pixel in map units) - F 20 GIS for Planetary Mappers

  21. Worldfile • Algebraic Form (six parameter affine transformation) x’ = Ax + By + C y’ = Dx + Ey + F where x’ = calculated x-coordinate of the pixel on the map y’ = calculated y-coordinate of the pixel on the map x = column number of a pixel in the image y = row number of a pixel in the image A = x-scale; dimension of a pixel in map units in x direction B,D = rotation terms (assumed to be zero) C,F = translation terms; x,y map coordinates of the center of the upper-left pixel E = negative of y-scale; dimension of a pixel in map units in y direction 21 GIS for Planetary Mappers

  22. Hands-on (Lambert Albedo – Exercise 05) Name: TES Albedo Filename: TES_Lambert_Albedo.png (Simple 8bit PNG format) Resolution: 8ppd Scale: 7.5kmpp Projection: Simple cylindrical, -180E to 180E, 90N to -90N, 'ocentric Layout: Single file Total Size: 2880x1440 pixels Details: Surface albedo. MGS/TES. 8 ppd/7.5km. Citation: Christensen et al., The Mars Global Surveyor Thermal Emission Spectrometer experiment: Investigation description and surface science results, J. Geophys. Res., 106, 23,823-23,871, 2001. from: http://www.mars.asu.edu/data/tes_albedo/ 22 GIS for Planetary Mappers

  23. Defining a GIS Worldfile Final worldfile looks like Worldfile for Lambert Albedo - file name = "TES_Lambert_Albedo.pgw" This PNG image is 2880 samples by 1440 lines and is from -180 to 180 longitude and -90 to 90 latitude (global) worldfile "TES_Lambert_Albedo.pgw" with descriptions: 0.125 // Xcellsize in degrees, 360 / num samples = 360 / 2880 0.0 // almost always 0 0.0 // almost always 0 -0.125 // Ycellsize usually = -X for square pixels -179.9375 // Upper left pixel (center) in X; -180 + (cellsize / 2) 89.9375 // Upper left pixel (center) in Y; 90 - (cellsize / 2) 0.125 0.0 0.0 -0.125 -179.9375 89.9375 23 GIS for Planetary Mappers

  24. PDS Worldfile • PDS uses same – but X,Y are in “pixel” space Worldfile (MOLA 4ppd megt90n000cb.lbl) OBJECT = IMAGE_MAP_PROJECTION ^DATA_SET_MAP_PROJECTION = "DSMAP.CAT" MAP_PROJECTION_TYPE = "SIMPLE CYLINDRICAL" A_AXIS_RADIUS = 3396.0 <KM> B_AXIS_RADIUS = 3396.0 <KM> C_AXIS_RADIUS = 3396.0 <KM> FIRST_STANDARD_PARALLEL = "N/A" SECOND_STANDARD_PARALLEL = "N/A" POSITIVE_LONGITUDE_DIRECTION = "EAST" CENTER_LATITUDE = 0.0 <DEGREE> CENTER_LONGITUDE = 180.0 <DEGREE> REFERENCE_LATITUDE = "N/A" REFERENCE_LONGITUDE = "N/A" LINE_FIRST_PIXEL = 1 LINE_LAST_PIXEL = 720 SAMPLE_FIRST_PIXEL = 1 SAMPLE_LAST_PIXEL = 1440 MAP_PROJECTION_ROTATION = 0.0 MAP_RESOLUTION = 4.0 <PIXEL/DEGREE> MAP_SCALE = 14.818 <KM/PIXEL> MAXIMUM_LATITUDE = 90.0 <DEGREE> MINIMUM_LATITUDE = -90.0 <DEGREE> WESTERNMOST_LONGITUDE = 0.0 <DEGREE> EASTERNMOST_LONGITUDE = 360.0 <DEGREE> LINE_PROJECTION_OFFSET = 360.5 SAMPLE_PROJECTION_OFFSET = 720.5 COORDINATE_SYSTEM_TYPE = "BODY-FIXED ROTATING" COORDINATE_SYSTEM_NAME = "PLANETOCENTRIC" END_OBJECT = IMAGE_MAP_PROJECTION 14818.0 (meters) 0.0 0.0 -14818.0 -10676369.0 X = SAMPLE_PROJ_OFFSET * MAP_SCALE * -1 5341889.0 Y = LINE_PROJ_OFFSET * MAP_SCALE http://pds-geosciences.wustl.edu/missions/mgs/megdr.html 24 GIS for Planetary Mappers

  25. Converting RAW FILES Example header (*.hdr) • ArcMap/GDAL • Create ESRI detached header • 8,16 bit, use extension *.bil or *.bsq • 32 bit file, use extension *.flt • Image & header must share filename • ArcMap • Can also use ERDAS detached header • GDAL • Can also use PCI Geomatic detached header ESRI Help: http://bit.ly/r3kIPJ GDAL Help: http://www.gdal.org/frmt_various.html#EHdr http://downloads.esri.com/support/whitepapers/other_/eximgav.pdf NCOLS xxx NROWS xxx XULCORNER xxx YULCORNER xxx CELLSIZE xxx NBITS 32 NODATA_VALUE xxx BYTEORDER <MSBFIRST | LSBFIRST>

  26. ASCII FILES • Regularly spaced • Add header to stream of “Z”s (filename *.asc) ESRI Help: http://bit.ly/r2GNFA NCOLS 480 NROWS 450 XULCORNER 378922 (or XLL) YULCORNER 4072345 CELLSIZE 30 NODATA_VALUE -32768 43 2 45 7 3 56 2 5 23 65 34 6 32 54 57 3 2 7 45 23 5 ...

  27. ASCII FILES • Irregularly spaced (randomly spaced) • Usually from a table (e.g. Lon, Lat, Value) • TAB or Comma delimited supported (*.tab, *.csv) • Once loaded (see next slide), you can then choose one of 10 interpolation methods in ArcMap • 3D Analyst, Spatial Analyst • GeoStatistical Analyst (interactive interpolation)

  28. ASCII FILES • Irregularly spaced (must create points prior to interpolation) http://bit.ly/psElPz

  29. ASCII FILES • Irregularly spaced (GMT, GDAL, QGIS, etc) 1.a) GMT Spherical interpolation http://gmt.soest.hawaii.edu/gmt/html/man/sphinterpolate.html# BlockMean or xyz2grd # http://www.soest.hawaii.edu/gmt/gmt/html/man/blockmean.html #set R=`minmax -I2 ascii.xyz` # Calculate the extent of the points #blockmeanascii.xyz -I0.01 -bo $R > temp.bm#If known extent set -Rxmin/xmax/ymin/ymax blockmean vesta_llr.txt -I0.0625 -bo -R0/360/-90/90 -:  > temp.bm # where -I resolution, use “-:” for lat,lonorder (leave off for lon, lat order)# where -bo means binary output and -bi means binary input (optional but faster)## run spherical interpolation (optionally run spherical TIN using sphtriangulate)# Spherical “Q1” =Smooth interpolation with local gradient estimates (more options avail.) sphinterpolatetemp.bm -Q1 -Gvesta_llr_sphInt_Q1.grd -I0.0625 -bo -R0/360/-90/90 -bi -:#now convert to GeoTiff or Raw (for import to ISIS using “raw2isis”)gdal_translate -of ENVI vesta_llr_sphInt_Q1.grd vesta_llr_sphInt_Q1.raw

  30. ASCII FILES • Irregularly spaced (GMT, GDAL, QGIS, etc) 1.b) GMT Cartesian interpolation (more typical - what MOLA//LOLA Team uses) http://www.soest.hawaii.edu/gmt/gmt/html/man/surface.html # BlockMean or xyz2grd # http://www.soest.hawaii.edu/gmt/gmt/html/man/blockmean.html #set R=`minmax -I2 ascii.xyz` # Calculate the extent of the points #blockmeanascii.xyz -I0.01 -bo $R > temp.bm #If known extent set -Rxmin/xmax/ymin/ymax blockmean vesta_llr.txt -I0.0625 -bo -R0/360/-90/90 -: > temp.bm # where -I resolution, use “-:” for lat,lon order (leave off for lon, lat order)# where -bo means binary output and -bi means binary input (optional but faster) #run spline interpolation (optionally run TIN using triangulate) surface temp.bm -Gvesta_llr_surface.grd -I0.0625 -bo -R0/360/-90/90 -bi -: #now convert to GeoTiff or Raw (for import into ISIS using raw2isis): gdal_translate -of ENVI vesta_llr_surface.grdvesta_llr_surface.raw

  31. Mosaic Raster Type A mosaic dataset is a collection of raster datasets (images) stored as a catalog & viewed as a dynamically mosaicked image.

  32. Mosaic data type • Demo Raster Riser

  33. Many and large image strategies Strategies for bringing in lots of images (e.g., ~100 projected NAC images @ 10-100 MB each) • Try Mosaics Data Type. • http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/What_is_a_mosaic_dataset/009t00000037000000/ • Related (see links at bottom): • See Raster Riser for instructions also (and interactive tool): • http://resources.arcgis.com/gallery/file/arcobjects-net-api/details?entryID=70A6CFB9-1422-2418-3473-F97D1590AFF1 • Also use caching basemaplayer and lock down scale options: • BaseMap Layers: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00s500000017000000.htm • Scale: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//006600000440000000.htm

  34. Virtual Image Functions Demo: Hillshade

  35. What is GDAL?(Geospatial Data Abstraction Library) • GDAL is a “translator library for raster geospatial data formats” • Open source • Used in many applications: GRASS, UMN MapServer, Google Earth, ArcGIS 9.x, etc. • Can handle many image formats for read and slightly fewer for write: AI Grid, Imagine, GeoTiff, JPEG, PNG, NetCDF, etc. From: Using Python, GDAL and NumPy for spatial analysis and modeling Outline

  36. GDAL (Geospatial Data Abstraction Library) • Presents an “abstract data model” for processing spatial data • Can be used directly from C/C++ and can be “wrapped” for use with Python, Perl, VB, C#, R, Java … • Early developers have chosen Python as their scripting language and documentation is relatively good for this. GDAL as of version 1.9.0 provides at least partial support for more than 120 raster geospatial data formats [ref] From: Using Python, GDAL and NumPy for spatial analysis and modeling Outline

  37. GDAL • Software that uses GDAL/OGR • World Wind Java NASA's open source virtual globe and world imaging technology • GRASS GIS • OSSIM • GvSIG • JMap • Quantum GIS • MapServer • Google Earth - A virtual globe and world imaging program. • OpenEV • SAGA GIS - a cross-platform open source GIS software • R - an open source statistical software with extensions for spatial data analysis • gdaltokmz, a Python module translating from GDAL-supported raster graphics formats to the Google Earth KMZ format • ArcGIS 9.2 can use GDAL for customized raster formats • TopoQuest - Internet topographic map viewer • Orfeo toolbox - A satellite image processing library • Biosphere3D – open source landscape scenery globe

  38. GDAL (Geospatial Data Abstraction Library) • Presents an “abstract data model” for processing spatial data • Can be used directly from C/C++ and can be “wrapped” for use with Python, Perl, VB, C#, R, Java … • Early developers have chosen Python as their scripting language and documentation is relatively good for this. From: Using Python, GDAL and NumPy for spatial analysis and modeling Outline

  39. GdalDEM– 8bit hillshades and Slope maps • Hillshade >gdaldemhillshade dem.img out_hillshade.tif -z 2       (z = exaggeration of 2) • colorize(using color.lut below) >gdaldem color-relief dem.img color.lut out_color.tif Merging colorized image and hillshade into a colorshade use: hsv_merge.py: http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ >hsv_merge.py out_color.tif out_hillshade.tifout_color-hillshade.tif For color mapping you need a mapping. Favorite so far is (rainbow). color.lut nv   0 0 0 // NoData 0%   purple 20%  blue 40%  aqua 60%  green 80%  yellow 100% red

  40. GdalDEM wrappers (@ ASTRO) • For LMMP color-defined shades/slopes run > gdal_colorshade_hsv.pl inDEM.cuboutClrShade.tif 2 • Where 2 is the exaggeration > gdal_slope_hsv.pl inDEM.cuboutSlope.tif2 Color legends will be copied to user’s directory but the ColorShade values will still need to be edited in Photoshop or other. Values will be written to screen.

  41. GdalDEM Results

  42. NumPy(Numerical Python) • An array/matrix package for Python • Well suited for image processing – i.e. one function can operate on the entire array • Slicing by dimensions and applying functions to these slices is concise and straightforward • Nearly 400 methods defined for use with NumPy arrays (e.g. type conversions, mathematical, logical, etc.) From: Using Python, GDAL and NumPy for spatial analysis and modeling Outline

  43. GDAL and NumPy • Since GDAL 1.3(?), GDAL has implemented NG (New Generation) Python bindings which includes NumPy • Process: Get raster band(s) Convert the raster band(s) to a NumPy array using ReadAsArray() Open GDALDataset Write out GDALDataset Process the raster band(s) using NumPy functionality Convert the NumPy array to GDAL raster bands using WriteAsArray() From: Using Python, GDAL and NumPy for spatial analysis and modeling Outline

  44. GDAL Virtual Format (*.vrt) • “The VRT a format driver for GDAL that allows a virtual dataset to be composed from other GDAL datasets withrepositioning, and algorithms potentially applied as well as various kinds of metadata altered or added. “ - Wikipedia

  45. Lazy raster processing with GDAL VRTs • Merge the tiles together • Reproject the merged DEM (using bilinear or cubic interpolation) • Generate the hillshade from the merged DEM • process-as-you-goimplementation: • >gdal_merge.py -of GTiff -o merged.tifdems*.tif • >gdalwarp-t_srs epsg:3310 -r bilinear -of GTiffmerged.tif merged_3310.tif • >gdaldemhillshade merged_3310.tif shade.tif -of Gtiff • lazy evaluation by using GDAL Virtual Rasters (VRT) • >gdalbuildvrtmerged.vrtdems*.tif • >gdalwarp-t_srs epsg:3310 -r bilinear -of VRT merged.vrt merged_3310.vrt • >gdaldemhillshade merged_3310.vrt shade.tif-of Gtiff • performs the intermediate steps, only creating GeoTiff as the final step http://blog.perrygeo.net/2010/02/18/lazy-raster-processing-with-gdal-vrts/

  46. Lazy raster processing with GDAL VRTs So what’s the advantage to doing it the VRT way? They both produce exactly the same output raster. Lets compare: http://blog.perrygeo.net/2010/02/18/lazy-raster-processing-with-gdal-vrts/

More Related