In [1]:
from datetime import datetime
from glob import glob
from IPython.display import display, HTML
from matplotlib import pyplot as plt
from osgeo import gdal, osr
from rasterio import Affine
from rasterio.rio import convert
from rasterio.warp import reproject, RESAMPLING
# local
from AlertaDengue.settings import (
MAPSERVER_URL, RASTER_PATH
)
from gis.geotiff import get_key_from_file_name
import fiona
import geopandas as gpd
import numpy as np
import os
import rasterio
import rasterio.plot
import rasterio.mask
# local
import AlertaDengue as alert_dengue
In [2]:
shp_path = '%s/static/shapefile' % os.path.dirname(alert_dengue.__path__[0])
raster_dir_path = RASTER_PATH
local_mapfile_dir = '/tmp/mapfiles/%s'
# mapserver variables
ms_shp_path = '/tiffs/'
ms_error_path = '/var/log/mapserver.log'
ms_cgi_path = MAPSERVER_URL + '?map=%s&'
ms_mapfile_name = '%s.map'
ms_mapfile_dir = '/maps/%s'
In [3]:
# Rio de Janeiro's conf
geocode = 3304557
city_name = 'Rio de Janeiro'
alert_level = 2
disease = 'dengue'
map_class = 'ndvi'
In [4]:
shapefile_path = os.path.join(shp_path, '%s.shp' % geocode)
gdf = gpd.GeoDataFrame.from_file(shapefile_path)
str_boundaries = [str(v) for v in gdf.bounds.iloc[0].values]
extent = ' '.join(str_boundaries)
extent_url = ','.join(str_boundaries)
crs_proj = gdf.crs['init']
wms_srs = crs_proj.upper()
In [5]:
data_range = {
'ndvi': (-2000.0, +10000.0),
'lst_day_1km': (0.0, 20000.0),
'lst_night_1km': (-30.0, 30.0),
'relative_humidity_2m_above_ground': (0.0, 100.0),
'specific_humidity_2m_above_ground': (0.0, 1.0),
'precipitation': (0, 200.0)
}
raster_classes = data_range.keys()
In [6]:
# mapfile templates
mapfile_template = '''
MAP
# The geographic extent (the rectangular area covered by the map) is
# defined by the keyword EXTENT. The rectangular area is specified by
# the coordinates of the opposite corners (the lower left and the upper
# right). These are coordinates of the southwest and the northeast corners
# EX:
# EXTENT -125.00 20.00 -65.00 50.00
# The geographic extent stretches from 125° west, 20° north to 65° west,
# 50° north.
CONFIG 'ON_MISSING_DATA' 'IGNORE'
CONFIG 'PROJ_LIB' './conf/'
CONFIG "MS_ERRORFILE" "%(ms_error_path)s"
CONFIG "CPL_DEBUG" "ON"
CONFIG "PROJ_DEBUG" "ON"
DEBUG 5
NAME "MAP_%(mapfile_name)s"
IMAGETYPE png
IMAGECOLOR 0 0 0
MAXSIZE 4000
SIZE 800 800
UNITS meters
EXTENT %(extent)s
OUTPUTFORMAT
NAME "png"
DRIVER AGG/PNG
MIMETYPE "image/png"
IMAGEMODE RGBA
EXTENSION "png"
FORMATOPTION "GAMMA=0.75"
END
OUTPUTFORMAT
NAME "GTiff"
DRIVER GDAL/GTiff
MIMETYPE "image/tiff"
IMAGEMODE RGBA
EXTENSION "tiff"
FORMATOPTION 'ATTACHMENT=%(mapfile_name)s.tiff'
END
OUTPUTFORMAT
NAME kml
DRIVER "KML"
MIMETYPE "application/vnd.google-earth.kml+xml"
IMAGEMODE RGBA
EXTENSION "kml"
FORMATOPTION 'ATTACHMENT=%(mapfile_name)s.kml'
END
PROJECTION
"init=%(crs_proj)s"
END
SHAPEPATH '%(shape_path)s'
WEB
METADATA
"wms_title" "Alerta Dengue"
"wms_onlineresource" "%(ms_cgi_path)s"
"wms_enable_request" "*"
"wms_srs" "%(wms_srs)s EPSG:3857"
"labelcache_map_edge_buffer" "-10"
# "wms_feature_info_mime_type" "text/html"
# "wms_format" "image/png"
END
IMAGEPATH '/tmp/map/'
IMAGEURL '/mapimg/'
END
%(include_layers)s
END
'''
mapfile_layer_template = '''
LAYER
NAME "%(city_name)s"
DATA "%(raster_path)s"
STATUS DEFAULT # ON
TYPE RASTER
OFFSITE 0 0 0
TEMPLATE "conf/template.html"
PROJECTION
"init=%(crs_proj)s"
END
METADATA
"wms_title" "%(city_name)s"
"wms_srs" "%(wms_srs)s EPSG:3857"
"wms_include_items" "all"
END
COMPOSITE
OPACITY 70
END # COMPOSITE
CLASS
NAME "class_hotcolors"
STYLE
RANGEITEM "style_hotcolors"
COLORRANGE 255 255 255 255 0 0
DATARANGE %(vmin)s %(vmax)s
END # STYLE
END # CLASS
END
'''
In [7]:
raster_file_path = os.path.join(
RASTER_PATH, 'meteorological', 'city', 'ndvi', '3304557', '20101219.tif'
)
shape_path = raster_dir_path
raster_name, src = '20101219.tif', rasterio.open(raster_file_path)
raster_file_name = '%s_%s' % (geocode, raster_name)
layer_name = '%s.map' % geocode
vmin, vmax = data_range[map_class]
layer_conf = {
'geocode': geocode,
'city_name': city_name.upper().replace(' ', '_'),
'layer_name': layer_name,
'rgb': '#FF9900',
'wms_srs': wms_srs,
'crs_proj': crs_proj,
'raster_path': raster_file_name,
'vmin': vmin,
'vmax': vmax,
'raster_title': map_class
}
In [8]:
# LAYER
include_layers = ''
layer_content = mapfile_layer_template % layer_conf
layer_relative_path = os.path.join(
'layers',
'meteorological',
map_class,
layer_conf['layer_name']
)
layer_path = local_mapfile_dir % layer_relative_path
include_layers += (
'INCLUDE "%s"' % layer_relative_path
)
layer_path_dir = os.path.dirname(layer_path)
In [9]:
# MAP
mapfile_name = '%s.map' % layer_conf['raster_title']
mapfile_path = local_mapfile_dir % mapfile_name
ms_mapfile_path = ms_mapfile_dir % mapfile_name
ms_config = {
'include_layers': include_layers,
'ms_error_path': ms_error_path,
'ms_cgi_path': ms_cgi_path % ms_mapfile_path,
'shp_path': ms_shp_path,
'extent': '-44.8893205505 -23.3689319629 -40.9585185182 -20.763205462',
'mapfile_name': 'NDVI',
'wms_srs': wms_srs,
'crs_proj': crs_proj,
'shape_path': ms_shp_path
}
mapfile_content = mapfile_template % ms_config
In [10]:
print('Generating ', mapfile_path, 'file ...')
print(mapfile_content)
In [11]:
print('Generating ', layer_path, 'file ...')
print(layer_content)