Skip to content

Instantly share code, notes, and snippets.

@epifanio

epifanio/ntg.map

Created March 14, 2025 18:28
Show Gist options
  • Save epifanio/88fd37162a4f584ff4e6cf3326545b0f to your computer and use it in GitHub Desktop.
Save epifanio/88fd37162a4f584ff4e6cf3326545b0f to your computer and use it in GitHub Desktop.
Download ZIP
Raw
ntg.map
MAP
NAME "Example MapServer WCS Multipart Issue"
CONFIG "MS_ERRORFILE" "stderr"
IMAGETYPE png
# EXTENT 190365.000 444075.000 190465.000 444175.000
# UNITS meters
# MAXSIZE 400
STATUS ON
DEBUG 5
WEB
METADATA
"ows_enable_request" "*"
"ows_fees" "NONE"
"ows_contactorganization" "Unknown"
"ows_schemas_location" "http://schemas.opengis.net"
# "ows_service_onlineresource" "https://wms.wps.met.no/get_wcs/test/wcs"
"ows_contactperson" "ContactCenter Unknown"
"ows_contactposition" "pointOfContact"
"ows_contactvoicetelephone" ""
"ows_contactfacsimiletelephone" ""
"ows_addresstype" ""
"ows_address" ""
"ows_city" "City"
"ows_stateorprovince" ""
"ows_postcode" ""
"ows_country" "Country"
"ows_contactelectronicmailaddress" "[email protected]"
"ows_hoursofservice" ""
"ows_contactinstructions" ""
"ows_role" ""
"ows_srs" "EPSG:4326 CRS:84"
"ows_accessconstraints" "otherRestrictions;http://creativecommons.org/publicdomain/mark/1.0"
END
END
# RESOLUTION 96
# DEFRESOLUTION 96
PROJECTION
"init=epsg:4326"
END
OUTPUTFORMAT
NAME "GEOTIFF_FLOAT32"
DRIVER "GDAL/GTiff"
MIMETYPE "image/tiff"
IMAGEMODE FLOAT32
EXTENSION "tif"
FORMATOPTION "NULLVALUE=-32768"
#FORMATOPTION "COMPRESS=DEFLATE"
#FORMATOPTION "PREDICTOR=3"
#FORMATOPTION "RESAMPLING=bilinear"
END # OUTPUTFORMAT
WEB
IMAGEPATH "/app/data/"
IMAGEURL "/app/data/"
METADATA
"ows_title" "Example MapServer WCS Multipart Issue"
"wcs_label" "Example MapServer WCS Multipart Issue"
"wcs_name" "Example MapServer WCS Multipart Issue"
"ows_abstract" "Abstract"
"wcs_description" "Abstract"
"wcs_fees" "NONE"
"wcs_keywordlist" ""
"wcs_metadatalink_type" "TC211"
"wcs_metadatalink_format" "application/vnd.ogc.csw.GetRecordByIdResponse_xml"
"wcs_metadatalink_href" "http://www.example.com?service=CSW&request=GetRecordById&version=2.0.2&outputSchema=http://www.isotc211.org/2005/gmd&elementSetName=full&ID=2d5d70b2-5b83-49c9-bf41-08413657a670"
"wcs_onlineresource" "https://wms.wps.met.no/get_wcs/test/wcs"
"ows_inspire_capabilities" "url"
"ows_languages" "dut"
"ows_inspire_dsid_code" "d3a939fa-2d2b-4eeb-8b23-14a5113b83eb"
"ows_inspire_dsid_ns" "http://example.com/"
"ows_inspire_metadataurl_format" "application/vnd.ogc.csw.GetRecordByIdResponse_xml"
"ows_inspire_metadataurl_href" "http:/www.example.com/geonetwork/srv/dut/csw-inspire?service=CSW&request=GetRecordById&version=2.0.2&outputSchema=http://www.isotc211.org/2005/gmd&elementSetName=full&ID=2d5d70b2-5b83-49c9-bf41-08413657a670"
END # METADATA
END # WEB
LAYER
NAME coverage
METADATA
"wcs_name" "coverage"
"wcs_label" "Title" # label describecoverage/1.0.0
"wcs_description" "Title" # title describecoverage/1.1.0, description describecoverage/1.0.0 - N.B. description en label zijn gelijk in describecoverage/1.1.0 lijkt bug in mapserver te zijn
"wcs_abstract" "Abstract" # abstract describecoverage/1.1.0
"wcs_srs" "EPSG:4326"
# "wcs_extent" "0.8715757 13.4142128 57.5062758 63.8639535"
# "wcs_resolution" "0.5 0.5"
"wcs_rangeset_name" "Rangeset Name" ### required to support DescribeCoverage request
"wcs_rangeset_label" "rangeset_label" ### required to support DescribeCoverage request
"wcs_formats" "GEOTIFF_FLOAT32"
"wcs_imagemode" "FLOAT32"
"wcs_bandcount" "1"
"wcs_band_names" "band_name"
"wcs_significant_figures" "1"
"hoogte_band_uom" "m"
"hoogte_band_definition" "http://example.com/phenomenon/definition"
"hoogte_band_description" "Band description"
"hoogte_interval" "0 200"
END # METADATA
TYPE RASTER ### required
STATUS ON
DATA /app/data/norway_nitrogen.tif
PROCESSING "RESAMPLE=BILINEAR"
PROJECTION
"init=epsg:4326"
END # PROJECTION
END # LAYER
END # MAP
Raw
wcs.py
from osgeo import gdal
import redis
from fastapi import Request, APIRouter, Query, HTTPException
from fastapi.responses import HTMLResponse, FileResponse, Response
import os
import sys
import logging
import logging.config
import yaml
import io
import mapscript
import xml.dom.minidom
import json
import multiprocessing
from urllib.parse import urlparse, parse_qs
import PIL.Image as Image
import uuid
sys.path.append('/app')
from pathlib import Path
import tempfile
import rasterio
from rasterio.mask import mask
from shapely.geometry import box
from fiona.crs import from_epsg
import pycrs
from rasterio.plot import show
import geopandas as gpd
import io
router = APIRouter()
# Load logging configuration from YAML file
with open('/app/logging.yaml', 'r') as file:
config = yaml.safe_load(file.read())
logging.config.dictConfig(config)
# Create logger
logger = logging.getLogger()
redis_client = redis.StrictRedis(host='redis', port=6379, db=0) # Create a Redis client
def bbox_clip(fp_in, minx, miny, maxx, maxy):
#fp_in = 'norway_nitrogen.tif'
#fp_out = 'norway_nitrogen_clipped.tif'
print('got here 01')
data = rasterio.open(fp_in)
print('got here 02')
# minx, miny, maxx, maxy = 9.073216441, 61.702432037999976, 10.288213072000019, 61.095326417
print(minx, miny, maxx, maxy)
bbox = box(minx, miny, maxx, maxy)
print('got here 03')
geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
print('got here 04')
def getFeatures(gdf):
"""Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
import json
return [json.loads(gdf.to_json())['features'][0]['geometry']]
coords = getFeatures(geo)
print('got here 05')
out_img, out_transform = mask(dataset=data, shapes=coords, crop=True)
print('got here 06')
out_meta = data.meta.copy()
print('got here 07')
# print(data.crs.data)
# epsg_code = int(data.crs.data['init'][5:])
try:
epsg_code = int(data.crs.data['init'][5:])
except Exception as e:
print(f"detect epsg code failed: {e}")
epsg_code = 4326
print('got here 08')
out_meta.update({"driver": "GTiff",
"height": out_img.shape[1],
"width": out_img.shape[2],
"transform": out_transform,
"crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4()})
print('got here 09')
#with rasterio.open(fp_out, "w", **out_meta) as dest:
# dest.write(out_img)
#with rasterio.open(io.BytesIO(), "w", **out_meta) as dest:
# dest.write(out_img)
# print(type(dest))
return (out_meta, out_img)
@router.get('/get_wcs/{data_id}/wcs', response_class=Response)
async def get_wcs(data_id: str, full_request: Request):
logging.debug(f"REQUEST URL: \n {full_request.url} \n")
map_object = mapscript.mapObj('/app/data/ntg.map')
layer = dict(full_request.query_params).get("LAYER")
layers = dict(full_request.query_params).get("LAYERS")
netloc = os.environ.get("NETLOC_ADDRESS", full_request.url.netloc)
scheme = os.environ.get("NETLOC_SCHEME", full_request.url.scheme)
ows_req = mapscript.OWSRequest()
ows_req.type = mapscript.MS_GET_REQUEST
try:
ows_req.loadParamsFromURL(str(full_request.query_params))
except mapscript.MapServerError:
logging.debug("loadParamsFromURL failed")
ows_req = mapscript.OWSRequest()
ows_req.type = mapscript.MS_GET_REQUEST
pass
if not full_request.query_params or ows_req.NumParams == 1:
logging.debug("Query params are empty. Force getcapabilities")
ows_req.setParameter("SERVICE", "WCS")
ows_req.setParameter("VERSION", "1.0.0")
ows_req.setParameter("REQUEST", "GetCapabilities")
else:
logging.debug(f"OWS REQUEST TYPE: {ows_req.type}")
logging.debug("Query params:")
logging.debug(f"NumParams: {ows_req.NumParams}")
for i in range(ows_req.NumParams):
logging.debug(f"{ows_req.getName(i)}: {ows_req.getValue(i)}")
if ows_req.getValueByName('REQUEST') not in ['GetCapabilities', 'DescribeCoverage', 'GetCoverage']:
logging.debug(f"got request: {ows_req.getValueByName('REQUEST')}")
dom = xml.dom.minidom.parseString(f"<p>{ows_req.getValueByName('REQUEST')}</p>")
result = dom.toprettyxml(indent="", newl="")
logging.debug(f"\n")
return Response(result, media_type='text/xml')
elif ows_req.getValueByName('REQUEST') == 'DescribeCoverage':
mapscript.msIO_installStdoutToBuffer()
dispatch_status = map_object.OWSDispatch(ows_req)
if dispatch_status != mapscript.MS_SUCCESS:
logging.debug(f"DISPATCH status: {dispatch_status}")
content_type = mapscript.msIO_stripStdoutBufferContentType()
mapscript.msIO_stripStdoutBufferContentHeaders()
_result = mapscript.msIO_getStdoutBufferBytes()
if content_type == 'application/vnd.ogc.wms_xml; charset=UTF-8':
content_type = 'text/xml'
logging.debug(f"content_type: {content_type}")
dom = xml.dom.minidom.parseString(_result)
result = dom.toprettyxml(indent="", newl="")
logging.debug(f"\n")
return Response(result, media_type=content_type)
elif ows_req.getValueByName('REQUEST') == 'GetCoverage':
format =ows_req.getValueByName('FORMAT')
print(format)
#mapscript.msIO_installStdoutToBuffer()
#dispatch_status = map_object.OWSDispatch(ows_req)
#result = mapscript.msIO_getStdoutBufferBytes()
#content_type="image/tiff"
#print("RESULTS Type: ", type(result))
# mapscript.msIO_stripStdoutBufferContentType()
#print(content_type)
#image = Image.open(io.BytesIO(result))
#image.save('/app/data/result.tif')
layer_name = ows_req.getValueByName('COVERAGE')
print(layer_name)
bbox_string = ows_req.getValueByName('BBOX')
print(bbox_string)
bbox = [float(i) for i in bbox_string.split(',')]
out_meta, out_img = bbox_clip('/app/data/norway_nitrogen.tif', bbox[0], bbox[1], bbox[2], bbox[3])
print(type(out_img))
print(type(out_meta))
# filename = uuid.uuid4()
tmp = tempfile.NamedTemporaryFile()
print('got here 10')
print(tmp.name)
with rasterio.open(tmp.name+'.tif', "w", **out_meta) as dest:
dest.write(out_img)
print(type(dest))
print('got here 11')
return FileResponse(tmp.name+'.tif', media_type="image/tiff")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment