Mapserver-Meteorological-Images


Mapserver - Meteorological Images


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

Setting variables


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()

Creating Mapfile for raster image


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)


Generating  /tmp/mapfiles/ndvi.map file ...

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" "/var/log/mapserver.log"
    CONFIG      "CPL_DEBUG" "ON"
    CONFIG      "PROJ_DEBUG" "ON"
    DEBUG       5
    
    NAME        "MAP_NDVI"
    
    IMAGETYPE   png
    IMAGECOLOR  0 0 0
    MAXSIZE     4000
    SIZE        800 800
    UNITS       meters
    EXTENT -44.8893205505 -23.3689319629 -40.9585185182 -20.763205462
       
    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=NDVI.tiff'
    END
    
    OUTPUTFORMAT
      NAME kml
      DRIVER "KML"
      MIMETYPE "application/vnd.google-earth.kml+xml"
      IMAGEMODE RGBA
      EXTENSION "kml"
      FORMATOPTION 'ATTACHMENT=NDVI.kml'
    END

    PROJECTION
        "init=epsg:4326"
    END
    
    SHAPEPATH '/tiffs/'

    WEB
      METADATA
        "wms_title" "Alerta Dengue"
        "wms_onlineresource" "http://172.17.0.2:80?map=/maps/ndvi.map&"
        "wms_enable_request" "*"
        "wms_srs" "EPSG:4326 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/meteorological/ndvi/3304557.map"
 
END


In [11]:
print('Generating ', layer_path, 'file ...')
print(layer_content)


Generating  /tmp/mapfiles/layers/meteorological/ndvi/3304557.map file ...

    LAYER
        NAME         "RIO_DE_JANEIRO"
        DATA         "3304557_20101219.tif"
        STATUS       DEFAULT # ON
        TYPE         RASTER
        OFFSITE      0 0 0
        
        TEMPLATE "conf/template.html"
        
        PROJECTION
          "init=epsg:4326"
        END
        
        METADATA
          "wms_title" "RIO_DE_JANEIRO"
          "wms_srs" "EPSG:4326 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 -2000.0 10000.0
         END # STYLE 

       END # CLASS 
    END