Skip to content

Jerrilyn's Work Mask Update

Christian Thomas edited this page Oct 16, 2019 · 5 revisions

Guide for Updating SkyTruth MTR Mask

FIPS Codes

Counties for the study area were selected using the following criteria, as provided by Appalachian Voices:

Are (1) within the "Central Appalachia" region as defined by the Appalachian Regional Commission

  • AND
  • (2a) Have reported surface coal production in at least one year since 1983
  • AND / OR
  • (2b) Have abandoned mine land sites in OSMRE's Abandoned Mine Land Inventory System (AMLIS)

See #55.

Here are the 74 counties’ FIPS codes.

Kentucky Tennessee Virginia West Virginia
21013 47001 51027 54005
21019 47013 51051 54011
21025 47025 51105 54015
21043 47035 51167 54019
21051 47049 51169 54025
21053 47129 51185 54039
21063 47133 51195 54043
21065 47137 51720 54045
21071 47141 54047
21089 47145 54053
21095 47151 54055
21109 54059
21115 54067
21119 54075
21121 54079
21125 54081
21127 54089
21129 54099
21131 54101
21133 54109
21135
21147
21153
21159
21165
21175
21189
21193
21195
21197
21199
21203
21205
21231
21235
21237

Hot Tip: for multiple counties in a row, select the first code then use shift-click to select the last code in the list. Also, watch out for Norton city: it is small.

Study Area

###Mask Data Acquisition

Use the TIGER/Line census data:https://www.census.gov/geo/maps-data/data/tiger-line.html

  • Download: AREAWATER, LINEARWATER, MIL, ROADS, UAC.
  • These are shapefiles of specific land uses/types

Hot Tip: Download Filezilla to speed up the download process. Using this application, you only need the url for the FTP download site for the data. Organize to use the Bash code below.

  • Unzip all of the downloads.

Hot Tip: Unzip the files en-mass using control-a then open with Archive Utility

Note: This site has land use data including a variety of uses not included in the previous downloads. We aren’t using this for the 2015 update, though.

These contain data of protected lands, like public forests and some military installments, and easement data, documents that prevent landowners from using land in certain ways

The study area is the (a lot of) counties, which you can select from the TIGER/line county file to create a reference shapefile. This is helpful for clipping the urban areas shapefile.

This shows the original study area, the darker region, overlayed with the updated region as of 2015.

###Merge

Hot Tip from Kevin: In Bash, use the script for SHP in tl*/.shp; do ogr2ogr -progress tl_2015_[STATE FIPS CODE][CLASS].shp “${SHP}” -append -nln tl_2015[STATE FIPS CODE]_[CLASS]; done to merge shapefiles. Modify based on file locations as needed. Substitute tl/*[identifier, eg areawater].shp depending on your file organization.

Remember, these must have the same geometries (all points or lines, etc).

In QGIS, check that the new shapefiles for roads and water (both kinds) match. Start by opening one file, then add the other layers: select the layer menu, add layer, add vector layer. Add the layer of military sites. Clip this to the study area layer. Repeat for urban areas file. Note that the one military site (as of 2015) does overlap with an urban area, so it is not included in the 2015 mask update.

###Buffer Save the layer and project to WGS84, Pseudo Mercator (3857). Merge files of similar geometries. Buffer each geometry 60 m. Merge buffers.

Hot Tip: Include the projection in the file name for clarity.

Hot Tip from Kevin: These buffers take a long time. Try this script: pip install fio #then fio buffer [file to be buffered] [file created] --distance 60 --buf-crs EPSG:3857 --jobs 4 [use all cores].

Here is the fully merged buffer for the 2015 study area:

###Rasterize Rasterize the merged buffer at 15m x 15m resolution.

Hot Tip: Once again, a very long process. Try: gdal_rasterize -burn 1 -l [full buffer file w/o extension] -a_nodata 0 -te -9558970.00000000000000 4250740.0000000000000000 -8622904.0000000000000000 4974270.0000000000000000 -tr 15 15 -ot Byte [input] [output] -a_srs EPSG:3857. Remember, the output is *.tif. This creates a raster with the same extents as the 2014 raster in the same projection, enabling us to compare results and locate differences.

(Run a Diff to locate the changes. Alternatively,) load in a google basemap (ie Google Satellite imagery) to check the validity of the new mask compared to the older mask. Some dirt roads are not featured in either the 2014 or ‘15 mask, so go through to get a better idea of what is and isn’t featured in each mask and what’s different. Double check with the old masks in the Earth Engine scripts.

Hot Tip: There shouldn’t be a lot of differences.

###Uploading to Earth Engine If you updated the study area, save as a keyhole markup language file by changing the file format (from .shp to .kml) within the save layer as window.

In Google Fusion Tables, upload the .kml version of the study area shapefile. In the Share window, update the settings so anyone with the link can view the table.

In Earth Engine, add the rasterized buffer as an asset through the options in the left side of the coding window.

Share the asset such that anyone can read the file. (This ensures the script incorporating this file runs for other users.)

Hot Tip from Kevin: This process takes a long time, try gdal_translate [raster file] [output] -co COMPRESS=DEFLATE -co TILED=YES -co INTERLEAVE=BAND -co NBITS=1 -co BIGTIFF=NO. Then add the new raster as an asset.

Update the code in the main window of the Earth Engine coding screen.

Hot Tip from Christian: If you updated the study area, add the new region using the format in line 5 and 7. Incorporate the new asset (raster) in line 13, add the new study area as a layer using the format in lines 23 and 24.

In the drop-down menu of the code window save option, choose to save as. Attribute to yourself, or don’t.

Upload the attributes to Google Drive, change the sharing settings to ‘Anyone with link can view’

###Troubleshooting

  1. Be careful applying bash scripts, use [command] --help for more information about the different flags within a command.
  2. “It’s the projection, stupid.” Relevant office poster. We ran into significant mismatching between the new and old raster not solved by reprojection. Be sure to save new files carefully, with new names that indicate the projection; stick to WGS 84, Pseudo mercator. Additionally, updating the shapefile of the study area is (occasionally) a pain. Check for sliver polygons when dissolving boundaries.
  3. With the rasterization and merging, the bash script is tricky. We produced multiple iterations that created files in World Mercator (54004) instead of 3857. Use spatialreference.org for help identifying other identifiers to try in the script.
  4. Make sure you share all components of the Earth Engine uploading so others can view and use your script updates.
  5. The blank raster image which all rasters should be based on is in 2015MTRmask file. Check that this aptly covers the western and southern bounds if used.