From d9fd71b094b34749cc771a7a2463421a74c99a25 Mon Sep 17 00:00:00 2001 From: Jordan A Caraballo-Vega Date: Tue, 11 Apr 2023 14:49:46 -0400 Subject: [PATCH] Adding folium_helper.py --- .gitignore | 3 + src/folium_helper.py | 129 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 132 insertions(+) create mode 100644 src/folium_helper.py diff --git a/.gitignore b/.gitignore index b6e4761..53f650f 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,9 @@ __pycache__/ *.py[cod] *$py.class +.DS_Store +*.hdf +*.gzip # C extensions *.so diff --git a/src/folium_helper.py b/src/folium_helper.py new file mode 100644 index 0000000..4a24a97 --- /dev/null +++ b/src/folium_helper.py @@ -0,0 +1,129 @@ +import folium +import glob +import numpy as np +import os +import rasterio as rio +import tempfile + +from rasterio.warp import calculate_default_transform, reproject, Resampling +from pyproj import Transformer + +# ----------------------------------------------------------------------------- +# Uses rasterio to open a raster, get the metadata and crs +# associated with it and get all the subdatasets in the file. +# This is very useful for hdf files such as MODIS hdfs. +# ----------------------------------------------------------------------------- +def print_subdatasets(filename): + bands_to_return = [] + with rio.open(filename) as dataset: + meta_data = dataset.meta + crs = dataset.read_crs() + + print([name for name in dataset.subdatasets if search_term in name]) + +# ----------------------------------------------------------------------------- +# Gets a tiff that has the correct metadata for that tile, gets the metadata +# from the source tif and copies to a destination tiff. +# ----------------------------------------------------------------------------- +def add_metadata_to_annual_product(filepath, model_type, year, tile): + metadata_pull_src = [fv for fv in glob.glob(os.path.join(filepath, "{}-1*-{}-MOD-*.tif".format(year, tile)))][0] + with rio.open(metadata_pull_src) as src: + src_meta = src.meta + dst_tiffs = [os.path.join(filepath, fn) for fn in os.listdir(filepath) if "{0}-{1}".format(year, tile) in os.path.basename(fn)] + [copy_meta(dst_tiff, src_meta, metadata_pull_src) for dst_tiff in dst_tiffs] + +# ----------------------------------------------------------------------------- +# Given a path to a tiff with no metadata, assign the metadata given to that +# tiff. +# ----------------------------------------------------------------------------- +def copy_meta(dst_path, src_meta, src_name): + print('Copying metadata from {} to {}'.format(src_name, dst_path)) + with rio.open(dst_path, 'r+') as dst: + dst.crs = src_meta['crs'] + dst.transform = src_meta['transform'] + +# ----------------------------------------------------------------------------- +# Given a tiff file as input, open the tiff and get the transform needed to +# reproject from the tiff's source crs to the one we want (EPSG:3857). +# For each band in the tiff, open then reproject it into the desired crs +# then write to a temporary file. Return the path to the temp file. +# ----------------------------------------------------------------------------- +def reproject_to_3857(input_tiff): + # Set desitnation CRS + dst_crs = f"EPSG:3857" + + # set out path + out_path_rproj = os.path.join(tempfile.gettempdir(), input_tiff.split('/')[-1].replace('.tif','-3857.tif')) + + with rio.open(input_tiff) as src: + # get src bounds and transform + transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds) + print('Transform: {}'.format(transform)) + print('Width: {} Height: {}'.format(width, height)) + kwargs = src.meta.copy() + kwargs.update({'crs': dst_crs, + 'transform': transform, + 'width': width, + 'height': height}) + + # reproject and write to file + with rio.open(out_path_rproj, 'w', **kwargs) as dst: + for i in range(1, src.count + 1): + reproject(source=rio.band(src, i), + destination=rio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=dst_crs, + resampling=Resampling.nearest) + return out_path_rproj + +# ----------------------------------------------------------------------------- +# In order for folium to work properly we need to pass it the bounding box +# of the tiff in the form of lat and lon. This is done by using rasterio. +# ----------------------------------------------------------------------------- +def get_bounds(tiff_3857): + with rio.open(tiff_3857) as src: + src_crs = src.crs['init'].upper() + min_lon, min_lat, max_lon, max_lat = src.bounds + bounds_orig = [[min_lat, min_lon], [max_lat, max_lon]] + bounds = [] + dst_crs = 'EPSG:4326' + for item in bounds_orig: + #converting to lat/lon + lat = item[0] + lon = item[1] + proj = Transformer.from_crs(int(src_crs.split(":")[1]), int(dst_crs.split(":")[1]), always_xy=True) + lon_n, lat_n = proj.transform(lon, lat) + bounds.append([lat_n, lon_n]) + center_lon = bounds[0][1] + (bounds[1][1] - bounds[0][1])/2 + center_lat = bounds[0][0] + (bounds[1][0] - bounds[0][0])/2 + return {'bounds': bounds, 'center': (center_lon, center_lat)} + +# ----------------------------------------------------------------------------- +# Use rasterio to open and read in the desired band name as a nd-array. +# ----------------------------------------------------------------------------- +def open_and_get_band(file_name, band_num=1): + with rio.open(file_name) as data: + b = data.read(band_num) + return b + +# ----------------------------------------------------------------------------- +# Given an nd-array (band) and the bounds in lat lon of the nd-array, return +# a folium layer. To add on the map. +# ----------------------------------------------------------------------------- +def get_overlay(band, meta_dict, name, opacity=1.0, show=True): + return folium.raster_layers.ImageOverlay(band, + bounds=meta_dict['bounds'], + name=name, + opacity=opacity, + show=show) + +# ----------------------------------------------------------------------------- +# We don't need to keep those temp files we made for the reprojections around. +# ----------------------------------------------------------------------------- +def cleanup(filename): + if os.path.exists(filename): + os.remove(filename) + else: + print('No file: {} exists.'.format(filename)) \ No newline at end of file