<jupyter_start><jupyter_text>GIS data processing  Shape files Load the required modules.<jupyter_code>from fiona import collection
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Create a fiona collection from the shape files.<jupyter_code>input = collection('Data/USA.gdb', 'r')<jupyter_output><empty_output><jupyter_text>Check the schema of the shape file.<jupyter_code>input.schema<jupyter_output><empty_output><jupyter_text>List the states in the US.<jupyter_code>states = set()
for region in input:
    states.add(region['properties']['NAME'])
print('number of states in US: {0}'.format(len(states)))
print('states: ', end='')
print(*states, sep=', ')<jupyter_output>number of states in US: 51
states: Maryland, Pennsylvania, Florida, Rhode Island, Mississippi, Hawaii, Arkansas, Massachusetts, South Dakota, Idaho, Colorado, New York, Arizona, Iowa, Michigan, Vermont, Wisconsin, Nebraska, California, New Hampshire, Oklahoma, Louisiana, Kentucky, Delaware, Connecticut, Alaska, Tennessee, Georgia, Indiana, Kansas, Oregon, South Carolina, Montana, Maine, Minnesota, Nevada, Ohio, Missouri, District of Columbia, New Mexico, North Dakota, Utah, Washington, Illinois, Alabama, North Carolina, Texas, West Virginia, New Jersey, Virginia, Wyoming<jupyter_text>Define a function that will return the points that define the polygon for a region.<jupyter_code>def get_region(input, n):
    return input[n]['geometry']['coordinates'][0][0]<jupyter_output><empty_output><jupyter_text>Plot all regions in the shapefile, this will be a map of the USA.<jupyter_code>plt.axes().set_aspect('equal')
for i in input.keys():
    region = get_region(input, i)
    x = [x[0] for x in region]
    y = [x[1] for x in region]
    plt.plot(x, y, 'b', linewidth=0.3)<jupyter_output><empty_output><jupyter_text>The bounds of the entire USA are given by an attribute.<jupyter_code>input.bounds<jupyter_output><empty_output><jupyter_text>Define a function that plots the state with the given name.<jupyter_code>def plot_state(input, name):
    plt.axes().set_aspect('equal')
    for i in input.keys():
        if input[i]['properties']['NAME'] == name:
            region = get_region(input, i)
            x = [t[0] for t in region]
            y = [t[1] for t in region]
            plt.plot(x, y, 'b', linewidth=0.7)<jupyter_output><empty_output><jupyter_text>Plot Florida's outline.<jupyter_code>plot_state(input, 'Florida')<jupyter_output><empty_output><jupyter_text>GeoJSON files  Reading GeoJSON files & Shapely Load the required modules.<jupyter_code>import geojson
import shapely
import shapely.geometry
data = geojson.load(open('Data/countries.geo.json'))
countries = data['features']
country_ids = [country['id'] for country in countries]
country_ids
def get_country_shape(countries, country_id):
    for country in countries:
        if country['id'] == country_id:
            return shapely.geometry.shape(country['geometry'])
usa = get_country_shape(countries, 'USA')
usa
canada = get_country_shape(countries, 'CAN')
canada<jupyter_output><empty_output><jupyter_text>GeoTIFF files Load the required modules.<jupyter_code>import gdal
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Open GeoTIFF file.<jupyter_code>data = gdal.Open('Data/cea.tif')<jupyter_output><empty_output><jupyter_text>Check the number of raster bands.<jupyter_code>data.RasterCount<jupyter_output><empty_output><jupyter_text>Check the number of pixels in the raster.<jupyter_code>data.RasterXSize, data.RasterYSize<jupyter_output><empty_output><jupyter_text>Retrieve the first (and only) band as a numpy array.<jupyter_code>raster_data = data.GetRasterBand(1).ReadAsArray()<jupyter_output><empty_output><jupyter_text>Note that the $X$ and $Y$ axes are swapped in the array.<jupyter_code>raster_data.shape<jupyter_output><empty_output><jupyter_text>Typically, the data in a band in a TIFF file are unsigned 8-bit integers.<jupyter_code>raster_data.dtype<jupyter_output><empty_output><jupyter_text>Plot the image, using a grayscale color map.<jupyter_code>plt.imshow(raster_data, cmap=cm.gray);<jupyter_output><empty_output><jupyter_text>Important features of a GeoTIFF file are its geo-transform, and its projection.<jupyter_code>data.GetGeoTransform()
data.GetProjection()<jupyter_output><empty_output><jupyter_text>To easily retrieve a GeoTIFF file's metadata, the `gdalinfo` command line utility is very useful.<jupyter_code>!gdalinfo Data/cea.tif
data.GetMetadata_Dict()<jupyter_output><empty_output><jupyter_text>To illustrate GeoTIFF file writing, we will normalize the raster band to the interval $\[0, 1\]$, and store it as single precisoin floating point data.<jupyter_code>new_raster_data = np.empty(raster_data.shape, dtype=np.float32)
new_raster_data = raster_data/raster_data.max()<jupyter_output><empty_output><jupyter_text>First, get the appropirate driver. GDAL supports many, but we need the `GTIFF` (i.e., GeoTIFF) driver.<jupyter_code>driver = gdal.GetDriverByName('GTIFF')<jupyter_output><empty_output><jupyter_text>Using the driver, a new file can be created. $X$ and $Y$ size, as well as the number of raster bands and the data type have to be specified.<jupyter_code>new_data = driver.Create('Data/gdal_to_be_removed.tif',
                         data.RasterXSize, data.RasterYSize, 1,
                         gdal.GDT_Float32, options=['COMPRESSION=LZW'])<jupyter_output><empty_output><jupyter_text>The projection and geo-transform are the same as for the original data.<jupyter_code>new_data.SetProjection(data.GetProjection())
new_data.SetGeoTransform(data.GetGeoTransform())<jupyter_output><empty_output><jupyter_text>Now we write the new numpy array to the first (and only) raster band.<jupyter_code>new_data.GetRasterBand(1).WriteArray(new_raster_data)<jupyter_output><empty_output><jupyter_text>Flush the data to disk.<jupyter_code>new_data.FlushCache()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11018780122984877, 'fraction_numerical': 0.003822502908426126, 'mean_word_length': 6.80544747081712, 'pattern_counts': {'":': 0, '<': 104, '<?xml version=': 0, '>': 104, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '427984', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2051, 'n_tokens_neox': 1995, 'n_words': 488, 'path': 'Python/Gis/gis.ipynb', 'repo_name': 'Gjacquenot/training-material', 'size': '6017'}