makepath / xarray-spatial Goto Github PK

• Flow Direction: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/how-flow-direction-works.htm
• Fill Depressions: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/creating-a-depressionless-dem.htm
• Flow Accumulation: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/how-flow-accumulation-works.htm
• Delineate watershed: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/watershed.htm

Part 1

Learn how to add arbitrary routes to a Panel applications.

• how do write a http handler in tornado?
• how do I connect that handler the running tornado instance in a Panel application (add question to panel github if needed https://github.com/holoviz/panel).

reference: https://docs.bokeh.org/en/latest/docs/dev_guide/server.html

Part 2

Implement the tile server using url:

• choice example dataset: Use data from the test data issue link
• implement custom route called /tile/{z}/{x}/{y}.png and returns Content-Type image/png.
• check if tile exists on disk
• if not, generate a single tile dynamically and save to disk with md5 of the url as the name (aefbc63827acb2991471098ab.png)
• add additional arguments for tile rendering.

Other

• add the panel zonal stats app example as a directory inside of examples examples/zonal-stats-panel
• add new example of using Tornado as a tile server.

The holy grail for this example is:

• super this command: zonal-stats -f dem.tif
  • This spawns a panel app which allows user to zoom and pan and view tiles of the file from the CLI command
  • Then user has drawing tools to create polygons
  • Polygons are sent to server and zonal stats are calculated off the original data at full resolution.

• Download validation data: http://xarray-spatial.s3.amazonaws.com/xarray-test-data.zip

For each one of the validations below:

1. Load the validation tif`
2. Calculate the equivalent output using Xarray-Spatial
3. Check if they are the same. If not, why not?

• Validate Hillshade against example output
• Validate Slope against example output
• Validate Aspect against example output
• Validate Curvature against example output
• Validate Zonal Statistics against example output
• Validate Viewshed against example output
• How are Edge pixels handled in validation data vs. Xarray-Spatial? Use slope as the example for your investigation.
• Validate Polygon Rasterization (datashader's) against the output in the test data

I think proximity with GREAT_CIRCLE looks still broken…can you try proximity with GREAT_CIRCLE using the example in user_guide.ipynb and confirm its good?

https://makepath.slack.com/files/UTLFF0HNH/F013YDJ2VUJ/screen_shot_2020-05-19_at_8.34.19_pm.png?origin_team=TTAAARNCR&origin_channel=CU8R1QL2V

https://earthenginepartners.appspot.com/science-2013-global-forest

Take Zonal Examples for User Guide, then expanded to have additional zonal functions

• zonal apply
• custom stats
• [ ]

https://openteams.com/projects/xarray-spatial

There is an interesting example of the approach implemented in Rust: https://github.com/mthh/contour-rs

Looks like the text and example don't match up. To fix these:

• get a local conda environment setup with xarray-spatial in editable mode git clone && cd && pip install -e .
• fix slope example
• fix curavture example

https://github.com/makepath/xarray-spatial/blob/master/xrspatial/zonal.py#L226

https://gis.apfo.usda.gov/arcgis/rest/services/NAIP/USDA_CONUS_PRIME/ImageServer

https://medium.com/geoai/high-resolution-land-cover-mapping-using-deep-learning-7126fee571dd

New Jupyter Notebook in examples called xarray-spatial-basic-ml.ipynb

The general template for ML examples:

• Use an Xarray-Spatial Analysis Function

• Take outputs and use as features for an scikit-learn model.

• Example 1 Clustering: Proximity To Pharmacies KMeans

• Regression Example

• Clustering Example using proximity and KMeans.

https://github.com/OSGeo/gdal/blob/master/gdal/alg/polygonize.cpp

• NDVI
• Move NDVI to new multispectral.py module
• List highest priorty multispectral tools
• Combine Bands
• ARVI (Atmospherically Resistant Vegetation Index)
• EVI (Enhanced Vegetation Index)
• EVI-2B (Enhanced Vegetation Index 2-band version)
• NBR (Normalized Burn Ratio)
• GCI (Green Chlorophyll Index)
• SIPI (Structure Insensitive Pigment Index)
• SAVI (Soil Adjusted Vegetation Index)

returns new xr.DataArray where we store the value of the nearest feature...instead of distance_to_nearest

from xrspatial.proximity import allocation
allocations = allocation(points_agg)
allocations.shape == points_agg.shape

We are going to port from PySAL and then give them appropriate credit...

• Quantile
• Equal-Interval
• Natural Breaks (1D K-means)

http://xarray.pydata.org/en/stable/related-projects.html

https://nrcs.app.box.com/v/elevation/folder/18389685693?page=3

• given a set a polygons, what resolution do I need to rasterize them to meet a minimum pixel count in each zone?
• we need to optimize for 4 things:
• xrange in data space - > The full x extent of the polygon GeoDataFrame
• yrange in data space - > The full x extent of the polygon GeoDataFrame

At what height, width do we are dropping polygons:

• width in pixel space
• height in pixel space

What is the smallest feature in the geodataframe by area?
If we were to represent that feature in a min of 25 pixels...what would the overall height width of the image be?

Code Idea

from xrspatial import suggest_zonal_canvas

my_polygon_df = GeoDataFrame(1M Polygons)
canvas_width, canvas_height = suggest_zonal_canvas(full_extent, min_feature_area, tolerance=.05)
cvs = ds.Canvas(width=canvas_width, height=canvas_height, ...)
cvs.polygons(my_polygon_df)

what should the size of my Canvas be to not drop polygons?

What should the canvas height width be so that the smallest polygon is rasterized with at least 25 pixels?

How about to not drop polygons and achieve a minimum pixel count?

Pharmacy Notebook has a decent example of this...

I want to be able to create the input dataset necessary for the ipyleaflet velocity layer.

Here is an example of what we want to get to:

from ipyleaflet import Map, TileLayer, basemaps
from ipyleaflet.velocity import Velocity
import xarray as xr
import os

if not os.path.exists('wind-global.nc'):
  url = 'https://github.com/benbovy/xvelmap/raw/master/notebooks/wind-global.nc'
  import requests
  r = requests.get(url)
  wind_data = r.content
  with open('wind-global.nc', 'wb') as f:
      f.write(wind_data)

center = [0, 0]
zoom = 1
m = Map(center=center, zoom=zoom, interpolation='nearest', basemap=basemaps.CartoDB.DarkMatter)

ds = xr.open_dataset('wind-global.nc')
display_options = {
    'velocityType': 'Global Wind',
    'displayPosition': 'bottomleft',
    'displayEmptyString': 'No wind data'
}

Let's consider adding https://pypi.org/project/tifffile/ as a dependency for tif io.

• add read_tif to readers.py which mirrors xr.open_rasterio

DEM data is stored as COG on S3

url = 'https://esip-pangeo-uswest2.s3-us-west-2.amazonaws.com/sciencebase/Southern_California_Topobathy_DEM_1m_cog.tif'
da = xr.open_rasterio(url)
da
<xarray.DataArray (band: 1, y: 225217, x: 316120)>
[71195598040 values with dtype=float32]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 3.825e+06 3.825e+06 3.825e+06 ... 3.6e+06 3.6e+06
  * x        (x) float64 1.795e+05 1.795e+05 1.795e+05 ... 4.956e+05 4.956e+05
Attributes:
    transform:   (1.0, 0.0, 179523.99999999822, 0.0, -1.0, 3824832.0)
    crs:         +init=epsg:26911
    res:         (1.0, 1.0)
    is_tiled:    1
    nodatavals:  (-3.4028230607370965e+38,)
    scales:      (1.0,)
    offsets:     (0.0,)
Check out the overviews:

src = rasterio.open(url, 'r')
[src.overviews(i) for i in src.indexes]
[[4, 8, 16, 32, 64, 128, 256, 512, 1023]]

taken from: https://nbviewer.jupyter.org/gist/rsignell-usgs/0f96bb9c0ca34a5dd0fc8131a7bbae1c

• Fork/Clone Xarray-Spatial code
• Create a new branch
• Create one 400px x 400px thumbnail for each analysis function using the crater lake DEM.
• For zonal tools, include geometry overlays in high-contrast color
• Add thumbnails to new images directory
• Update markdown below with image paths
• Add markdown to readme and create Pull Request

Example Markdown for Image Grid

DEM	SLOPE	CURVATURE
HILLSHADE	ZONAL

• create new directory /examples/user_guide/
• add /examples/user_guide/0_Getting_Setup.ipynb
• add /examples/user_guide/1_Surface.ipynb
• add /examples/user_guide/2_Proximity.ipynb
• add /examples/user_guide/3_Zonal.ipynb
• add /examples/user_guide/4_Focal.ipynb
• add /examples/user_guide/5_Classification.ipynb
• add /examples/user_guide/6_Conversion.ipynb ()
• add /examples/user_guide/7_Tiles.ipynb (Generate tileset)
• add /examples/user_guide/8_Remote_Sensing.ipynb (NDVI, Contrast Enhancment, Segmentation)

https://nmbbmapping.org/arcgis/rest/services/AllProviders/MapServer

• get cartoframes api key from Brendan
• create new notebook based on pharmacy notebook to demonstrate integration with cartoframes

• Let's test and move over code from datashader PR here:
  https://github.com/tomwhite/datashader/tree/tile-y-coordinate-fix
• Test with S3 Backend
• Test with World Polygons Examples
• Use span to set global shading min/ max range
• add wgs84 tile definition
• PEP8 compliance
• add to_ogc_tile_metadata which outputs a tile GetCapabilities.xml file: http://schemas.opengis.net/wmts/1.0/profiles/wmts-simple/examples/wmtsGetCapabilities_response_OSM.xml
• add to_esri_tile_metadata which outputs a tile MapServer file (json):

{
"serviceDescription" : "Test Map Service Description",
"mapName" : "Street Map Pro Data",
"description": "Street Map USA",
"copyrightText" : "ESRI",
"layers": [
    {"id" : 0, "name" : "Cities", "defaultVisibility" : true, "parentLayerId" : -1, "subLayerIds" : null},
    {"id" : 1, "name" : "States", "defaultVisibility" : true, "parentLayerId" : -1, "subLayerIds" : null},
    {"id" : 2, "name" : "Counties", "defaultVisibility" : false, "parentLayerId" : -1, "subLayerIds" : [3, 4]},
    {"id" : 3, "name" : "Large Counties", "defaultVisibility" : false, "parentLayerId" : 2, "subLayerIds" : null},
    {"id" : 4, "name" : "Small Counties", "defaultVisibility" : false, "parentLayerId" : 2, "subLayerIds" : null}
],
"spatialReference" : {"wkid" : 4326},
"singleFusedMapCache" : true,
"tileInfo": {
    "rows" : 512, "cols" : 512, "dpi" : 96, "format" : "JPEG", "compressionQuality" : 75,
    "origin" : {"x" : -130.0, "y" : 50.0},
    "spatialReference" : {"wkid" : 4326},
    "lods": [
        {"level" : 0, "resolution" : 8.46, "scale" : 32000.0 },
        {"level" : 1, "resolution" : 4.23, "scale" : 16000.0 },
        {"level" : 2, "resolution" : 2.11, "scale" : 8000.0 },
        {"level" : 3, "resolution" : 1.05, "scale" : 4000.0 },
        {"level" : 4, "resolution" : 0.52, "scale" : 2000.0 }
    ]
},
"initialExtent" : {
    "xmin" : -109.55, "ymin" : 25.76, "xmax" : -86.39, "ymax" : 49.94, 
    "spatialReference" : {"wkid" : 4326}
},
"fullExtent" : {
    "xmin" : -130.0, "ymin" : 24.0, "xmax" : -65.0, "ymax" : 50.0,
    "spatialReference" : {"wkid" : 4326}
},
"units" : "esriDecimalDegrees",
"supportedImageFormatTypes": "PNG32,PNG24,PNG,JPG,DIB,TIFF,EMF,PS,PDF,GIF,SVG,SVGZ",
"documentInfo": {
    "Title" : "StreetMap USA.mxd",
    "Author" : "ESRI Data Team",
    "Comments" : "ESRI Data and Maps 2004",
    "Subject" : "Street level data for the US",
    "Category" : "vector",
    "Keywords" : "StreetMap USA"
}
}```

- [ ] add mbtile output

• Heat Index

Support for dask arrays will allow for distributed spatial analysis.

Convolution filter dask array pseudo example

def distributed_hillshade(block):
  return hillshade(block)

g = da.overlap.overlap(mydaskarray, depth={0: 2, 1: 2},
                    boundary={0: 'periodic', 1: 'periodic'})
g2 = g.map_blocks(distributed_hillshade)
result = da.overlap.trim_internal(g2, {0: 2, 1: 2})

Started a repo here for it: https://github.com/makepath/xarray-spatial-feedstock

Right now proximity requires x y named fields in the xarray.DataArray. Lets expose x_field and y_field optional arguments as overrides.

Add new kernels:

• add custom (allow user to supply 2d array representing convolution kernel)
• add circle
• add wedge
• add cone

I am curious about the project and personally find that web accessible docs are very helpful

Looks like @perrygeo has an awesome library for this...we just want it implemented in numba...@thuydotm let's make sure to give attribution
https://github.com/perrygeo/jenks/blob/master/jenks.pyx

• update release notes
• add conda forge feedstock
• deploy package to pypi
• deploy package to anaconda.org

returns new xr.DataArray where we store the direction (0.0-259.0 degrees) to the nearest feature...instead of distance_to_nearest

from xrspatial.proximity import direction
directions = direction(points_agg)
directions.shape == points_agg.shape

Reproduction:

import datashader as ds
from xrspatial import proximity

df1 = pd.DataFrame(dict(x=[-74.39], y=[43.0], val=[122]))
canvas = ds.Canvas(plot_height=300, plot_width=300, x_range=(-74, -75), y_range=(42,44))
df1_agg = canvas.points(df1, 'x', 'y', ds.sum('val'))
df1_agg.data = df1_agg.data.astype('float64')
df1_agg.data[np.isnan(df1_agg.data)] = 0.0
df1_proximity = proximity(df1_agg, distance_metric='GREAT_CIRCLE')
np.nanmin(df1_proximity)

• add initial code
• add tests
• add docs
• add example
• add tutorial