# Data accessed using command:

wget https://naciscdn.org/naturalearth/10m/raster/NE2_HR_LC_SR_W.zip -O {output_dir}/NE2_HR_LC_SR_W.zip

# Data processed using the following steps:

unzip {input_dir}/*.zip -d {output_dir} 

gdalwarp -t_srs EPSG:3413 -r bilinear -tr 500 500 -te -5774572.727595 -5774572.727595 5774572.727595 5774572.727595 -dstnodata 0 -wo SOURCE_EXTRA=100 -wo SAMPLE_GRID=YES {input_dir}/NE2_HR_LC_SR_W.tif {output_dir}/warped.tif

gdalwarp -cutline {assets_dir}/latitude_shape_40_degrees.geojson -crop_to_cutline -co COMPRESS=DEFLATE {input_dir}/warped.tif {output_dir}/warped_and_cut.tif

cp {input_dir}/warped_and_cut.tif {output_dir}/input.tif && gdal_calc.py --calc="numpy.invert(numpy.isnan(A)).astype(int)" --outfile={output_dir}/mask.tif -A {output_dir}/input.tif --A_band=1

cp {input_dir}/* {output_dir}/ && gdal_translate -b 1 {input_dir}/input.tif {output_dir}/b1.tif && gdal_translate -b 2 {input_dir}/input.tif {output_dir}/b2.tif && gdal_translate -b 3 {input_dir}/input.tif {output_dir}/b3.tif

gdal_merge.py -o {output_dir}/merged.tif -co COMPRESS=DEFLATE -separate {input_dir}/b1.tif {input_dir}/b2.tif {input_dir}/b3.tif {input_dir}/mask.tif

gdal_translate -co TILED=YES -co COMPRESS=JPEG -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR -b 1 -b 2 -b 3 -mask 4 --config GDAL_TIFF_INTERNAL_MASK YES {input_dir}/merged.tif {output_dir}/compressed.tif

cp {input_dir}/compressed.tif {output_dir}/overviews.tif && gdaladdo -r average {output_dir}/overviews.tif 2 4 8 16