1

I am trying to establish a data pipeline that loop through a series of bounding boxes (defined from another dataset) and then looks for the intersection of a "fire boundary" with "populated areas".

I'm able to render the census data I care about as follows -- we only care about the yellow, orange, and red as defined....

dataset = ee.FeatureCollection('TIGER/2010/Blocks')
palette = ['black', 'yellow', 'orange', 'red']
opacity = 0.7
visParams = {
  'min': 30.0,
  'max': 100.0,
  'palette': palette,
  'opacity': opacity,
  'region': LAWarm_AOI}

# Turn the strings into numbers

popImage = ee.Image().float().paint(dataset, 'pop10');

my_map.addLayer(popImage, visParams, 'TIGER')

my_map

we're then looking to get the overlap with a fire boundary as visualized in pink here...

#confirm the burn area is in that vicinity!
MODIS = ee.ImageCollection("MODIS/006/MCD64A1").filter(ee.Filter.date(getDatePre(fire_perims,i)[0],getDatePost(fire_perims,i)[1]))
burnedArea = MODIS.select('BurnDate')

palette = ['pink']
burnedArea_parameters = {'min': 0,
                   'max': 1,
                   'dimensions': 512,
                   'palette': palette,
                   'region': LAWarm_AOI}
my_map.addLayer(burnedArea, burnedArea_parameters, 'burnyMBF')
my_map

I'm able to calculate an intersection using the following....

fireSingleImage=burnedArea.reduce(ee.Reducer.median())
fireVector=fireSingleImage.reduceToVectors()
popVector=popImage.reduceToVectors()
fireGeo=fireVector.geometry()
popGeo=popVector.geometry()
intersection = fireGeo.intersects(popGeo, ee.ErrorMargin(1))
intersection = ee.Geometry.Polygon(intersection)
print(intersection)

but when I go to visualize this intersection to confirm I'm picking up the right area as so...

Intersectpalette = ['green']
Intersect_parameters = {'min': 0,
                   'max': 1,
                   'dimensions': 512,
                   'palette': Intersectpalette,
                   'region': LAWarm_AOI}
my_map.addLayer(intersection, Intersect_parameters, 'intersection')
my_map

I get back an error that says EEException: Image.reduceToVectors: The default WGS84 projection is invalid for aggregations. Specify a scale or crs & crs_transform.

Again... this is for a data pipeline so I need to be able to execute this routine several times... ultimately we're just trying to get the intersection geometries for ~17,000 examples like above and store them for further processing

Any thoughts on how to establish a routine to do this?

FULL CODE BELOW:

#Requirement.... 

#Python 3.6 "Earth Engine" environment w. following installs.... 
#------- Jupyter Notebook
#------- Earth Engine
# Developer docs here: https://developers.google.com/earth-engine
import ee
ee.Authenticate() #should take you to Google SSO page --- use internal credentials!
ee.Initialize()

#------- Folium for interactive maps
import folium
#------- Geehydro to help w/ JS emulation for interactive maps
import geehydro
#------- Datetime just to help extract metadata
import datetime as dt
#------- Image from IPython to help w. displaying data
from IPython.display import Image
#pandas for data frames
import pandas as pd
#np for matrices, etc
import numpy as np
#parse NOAA datetimes
from datetime import datetime
#json compatability
import json
#web requests
import requests
#geopandas
import geopandas as gpd
#geolocations
from owslib.wms import WebMapService
#pyplot
import matplotlib.pyplot as plt
#shapely
import shapely   
#time
import time
#os
import os
#cv2 for OpenCV computer vision
import cv2    
#math for degree to rad
import math

#get everything into one nice geodataframe
fire00 = gpd.read_file('2000_perimeters_dd83.shp')
fire01 = gpd.read_file('2001_perimeters_dd83.shp')
fire02 = gpd.read_file('2002_perimeters_dd83.shp')
fire03 = gpd.read_file('2003_perimeters_dd83.shp')
fire04 = gpd.read_file('2004_perimeters_dd83.shp')
fire05 = gpd.read_file('2005_perimeters_dd83.shp')
fire06 = gpd.read_file('2006_perimeters_dd83.shp')
fire07 = gpd.read_file('2007_perimeters_dd83.shp')
fire08 = gpd.read_file('2008_perimeters_dd83.shp')
fire09 = gpd.read_file('2009_perimeters_dd83.shp')
fire10 = gpd.read_file('2010_perimeters_dd83.shp')
fire11 = gpd.read_file('2011_perimeters_dd83.shp')
fire12 = gpd.read_file('2012_perimeters_dd83.shp')
fire13 = gpd.read_file('2013_perimeters_dd83.shp')
fire14 = gpd.read_file('2014_perimeters_dd83.shp')
fire15 = gpd.read_file('2015_perimeters_dd83.shp')
fire16 = gpd.read_file('2016_perimeters_dd83.shp')
fire17 = gpd.read_file('2017_perimeters_dd83.shp')
fire18 = gpd.read_file('2018_perimeters_dd83.shp')
frames = [fire18, fire17, fire16, fire15, fire14, fire13, fire12, fire11, fire10, fire09, fire08, fire07, fire06, fire05, fire04, fire03, fire02, fire01, fire00]
all_fires = gpd.GeoDataFrame( pd.concat( frames, ignore_index=True) )
len(all_fires)

#trim to just California fires... note that below code actually call all fires 
CAall_fires = gpd.GeoDataFrame(all_fires[all_fires.state=='CA'])
len(CAall_fires)

#visualize just to get an idea of the shape of the data... "uniquefire", "perimeterdate" and "geometry" are key fields
fire_perims = CAall_fires.copy()
data_top = fire_perims.head()  
data_top

#Helper functions for manipulating degrees to radians, radians to degrees etc... 

# degrees to radians
def deg2rad(degrees):
    return math.pi*degrees/180.0
# radians to degrees
def rad2deg(radians):
    return 180.0*radians/math.pi

# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
def WGS84EarthRadius(lat):
    # http://en.wikipedia.org/wiki/Earth_radius
    An = WGS84_a*WGS84_a * math.cos(lat)
    Bn = WGS84_b*WGS84_b * math.sin(lat)
    Ad = WGS84_a * math.cos(lat)
    Bd = WGS84_b * math.sin(lat)
    return math.sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )

# Semi-axes of WGS-84 geoidal reference
WGS84_a = 6378137.0  # Major semiaxis [m]
WGS84_b = 6356752.3  # Minor semiaxis [m]

def boundingBox(longitudeInDegrees, latitudeInDegrees, halfSideInKm):
    lon = deg2rad(longitudeInDegrees)
    lat = deg2rad(latitudeInDegrees)
    halfSide = 1000*halfSideInKm

    # Radius of Earth at given latitude
    radius = WGS84EarthRadius(lat)
    # Radius of the parallel at given latitude
    pradius = radius*math.cos(lat)

    min_lat = lat - halfSide/radius
    max_lat = lat + halfSide/radius
    min_lon = lon - halfSide/pradius
    max_lon = lon + halfSide/pradius
    Bbox = (rad2deg(min_lon), rad2deg(min_lat), rad2deg(max_lon), rad2deg(max_lat))
    return Bbox


#
##
###
#### manually selected i=1051, Camp Fire (Paradise) from 2018 for testing purposes
###
##
#
#
##
###
#### ToDo: May need to dedupe fire boundaries... for example, there are multiple Camp fires (use "incidentname")
###
##
#
#testing the bounding box approach (fire should be in Paradise Area).... will later loop through i.... 
i = 1051
distance = 2

def getbox(i, distance): 
    if type(all_fires['geometry'][i]) == shapely.geometry.polygon.Polygon:
        test_polygon = all_fires['geometry'][i]
        ade_center = np.array(test_polygon.representative_point()) 
        polybbox = boundingBox(ade_center[0], ade_center[1], distance)
    elif type(all_fires['geometry'][i]) == shapely.geometry.multipolygon.MultiPolygon:
        test_polygon = all_fires['geometry'][i][0]
        ade_center = np.array(test_polygon.representative_point()) 
        polybbox = boundingBox(ade_center[0], ade_center[1], distance)
    return polybbox

getbox(i,distance)

#make sure you're getting good lat long
print(getbox(i,distance))

#pull up a map to check is centered over Paradise California Area
myade = np.array(all_fires['geometry'][i].representative_point())
my_map=folium.Map(location=[myade[1], myade[0]], zoom_start=10)
my_map

#these helper functions will give you a series of dates related to each fire as you loop through (note i above)
#dates help to make sure you're getting right fire boundaries and relevant before/after imagery

fire_perims = all_fires.copy()

#find the dates for the after photo --- currently set to 20 days later... 
def getDatePost(fire_perims, i):
    #get date and turn it into date time... 
    date = fire_perims['perimeterd'][i]
    firedate = dt.datetime.strptime(date, '%Y-%m-%d')
    firedate_str = firedate.strftime("%Y-%m-%d")
    interval_enddate = firedate + dt.timedelta(days = 20)
    interval_enddate_str = interval_enddate.strftime("%Y-%m-%d")
    return (firedate_str,interval_enddate_str)

#find the dates for the before photo... 
def getDatePre(fire_perims, i):
    #get date and turn it into date time... 
    date = fire_perims['perimeterd'][i]
    firedate = dt.datetime.strptime(date, '%Y-%m-%d')
    interval_startdate = firedate - dt.timedelta(days = 100)
    firedate_str=interval_startdate.strftime("%Y-%m-%d")
    interval_enddate= firedate - dt.timedelta(days = 20)
    interval_enddate_str=interval_enddate.strftime("%Y-%m-%d")
    return (firedate_str,interval_enddate_str)

print(fire_perims['perimeterd'][i])
print(getDatePost(fire_perims,i)[0])
print(getDatePost(fire_perims,i)[1])
print(getDatePre(fire_perims,i)[0])
print(getDatePre(fire_perims,i)[1])

#Initialize landsat as all Landsat 8 Surface Reflectance Tier 1 filtered for the three months preceding the fire
#NOTE: this is where Planet or Sentinel may come in handy!!! 
landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
landsat = landsat.filterDate(getDatePre(fire_perims,i)[0],getDatePre(fire_perims,i)[1]) 
#Set an AOI (Area of Interest)
LAWarmbox = getbox(i,distance)
LAWarm_AOI = ee.Geometry.Rectangle([ LAWarmbox[0], LAWarmbox[1],
                                       LAWarmbox[2], LAWarmbox[3]])

#call the GeoJSON
landsat_AOI = landsat.filterBounds(LAWarm_AOI)

#Check on the number of images and bands available in the time of interest
print('Total number:', landsat_AOI.size().getInfo())
landsat_AOI.first().bandNames().getInfo()

#This looks at the perimeter we've defined above and the dates we've defined above to make sure there is a fire there
#more of a helpful visual confirmation vs. something for data pipeline
#should be a pink area in the middle of the visual... 

MODIS = ee.ImageCollection("MODIS/006/MCD64A1").filter(ee.Filter.date(getDatePre(fire_perims,i)[0],getDatePost(fire_perims,i)[1]))
burnedArea = MODIS.select('BurnDate')

palette = ['pink']
burnedArea_parameters = {'min': 0,
                   'max': 1,
                   'dimensions': 512,
                   'palette': palette,
                   'region': LAWarm_AOI}
my_map.addLayer(burnedArea, burnedArea_parameters)
my_map

#get the least cloudy image and print out the cloud cover percentage
#alternatively could use cloud cover removal but this feels good enough for our purposes

least_cloudy = ee.Image(landsat_AOI.sort('CLOUD_COVER').first())
print('Cloud Cover (%):', least_cloudy.get('CLOUD_COVER').getInfo())

#
date = ee.Date(least_cloudy.get('system:time_start'))
time = date.getInfo()['value']/1000
dt.datetime.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')

#load least cloudy image from before the fire... may be a partial visual
sat_parameters = {'min': 100,
              'max': 1000,
              'dimensions': 1000,
              'bands': ['B4','B3', 'B2'],   #these are RGB bands!!!
              'region': LAWarm_AOI}
my_map.addLayer(least_cloudy, sat_parameters)
my_map

#TIGER data set goes and gets census tracks and makes any with >30 Yellow, Orange, and Red... 
#so we are targeting area where yellow, orange, and red overlap
#note: can take a few seconds to lod... opacity helps confirm overlap... 

dataset = ee.FeatureCollection('TIGER/2010/Blocks')
palette = ['black', 'yellow', 'orange', 'red']
opacity = 0.7
visParams = {
  'min': 30.0,
  'max': 100.0,
  'palette': palette,
  'opacity': opacity,
  'region': LAWarm_AOI}

# Turn the strings into numbers

popImage = ee.Image().float().paint(dataset, 'pop10');

my_map.addLayer(popImage, visParams, 'TIGER')

my_map

#
##
###
#### THIS IS WHERE INTERSECTION CODE STARTS... 
#### Note the error likely needs to be corrected using CRS ("Coordinate Refrence System") vvv
#### <Follow ME>. https://www.earthdatascience.org/workshops/gis-open-source-python/reproject-vector-data-in-python/
###
##
#
fireSingleImage=burnedArea.reduce(ee.Reducer.median())

#fireVector and popVector are fire boundary and tiger population layer respectively as Vectors... then geometryiies
fireVector=fireSingleImage.reduceToVectors()
popVector=popImage.reduceToVectors()
fireGeo=fireVector.geometry()
popGeo=popVector.geometry()

#intersection is then calculated and converted to a Polygon... 
intersection = fireGeo.intersects(popGeo, ee.ErrorMargin(1))
intersection = ee.Geometry.Polygon(intersection)
print(intersection.getInfo())

#should yield an error that says... "Image.reduceToVectors: The default WGS84 projection is invalid 
#for aggregations. Specify a scale or crs & crs_transform.

#In theory if intersection is properly defined... this will then make the right area light up green... 
Intersectpalette = ['green']
Intersect_parameters = {'min': 0,
                   'max': 1,
                   'dimensions': 512,
                   'palette': Intersectpalette,
                   'region': LAWarm_AOI}
my_map.addLayer(intersection, Intersect_parameters, 'intersection')
my_map
1

First, a tip to increase the chance someone will help you out with your questions. Before submitting, try to debug it as far as possible yourself. Narrow down where things go wrong. Remove all code not needed to reproduce your problem, but make sure the code still is executable for someone else, if at all possible.

There are a couple of problems with the code. Here's part of your script with inline comments that should sort it out:

import ee

ee.Initialize()

# Test region
region = ee.Geometry.Polygon([[
    [12.34699844239259, 42.01815395876468],
    [12.34699844239259, 41.64158937004622],
    [12.85236953614259, 41.64158937004622],
    [12.85236953614259, 42.01815395876468]
]])

dataset = ee.FeatureCollection('TIGER/2010/Blocks')
popImage = ee.Image().float().paint(dataset, 'pop10')

#confirm the burn area is in that vicinity!
MODIS = ee.ImageCollection("MODIS/006/MCD64A1") \
    .filterDate('2019-01-01', '2019-02-01')
burnedArea = MODIS.select('BurnDate')
fireSingleImage = burnedArea.reduce(ee.Reducer.median())
# Requirements for reduceToVectors():
# - Image must have an int data type. 
#   I haven't investigated what range of values you have in fireSingleImage and popImage. 
#   If they're small, you have to scale them up before casting to an int
# - Image must be bounded (clip it), or you have to provide a geometry
# - A scale, or crs and crsTransformation must be specified
fireVector = fireSingleImage.int().reduceToVectors(geometry=region, scale=463)
# Add some debug logs to find out exactly where and how things break. 
# Remember to call getInfo() on any EE objects.
print("fireVector = {}".format(fireVector.getInfo())) 
popVector = popImage.int().reduceToVectors(geometry=region, scale=463)
print("popVector = {}".format(popVector.getInfo()))
fireGeo = fireVector.geometry()
popGeo = popVector.geometry()
# intersects() returns a boolean - True if the two geometries intersects, otherwise False
# You want to use intersection() here instead
# intersection = fireGeo.intersects(popGeo, ee.ErrorMargin(1))
intersection = fireGeo.intersection(popGeo, ee.ErrorMargin(1))
# The intersection doesn't have to be a simple polygon, so it can fail when you try to cast it to one.
# It's already a geoemtry, and it doesn't really matter which type. Just use it
# intersection = ee.Geometry.Polygon(intersection) 
print('intersection = {}'.format(intersection.getInfo()))
2
  • thank @Daniel ... the full code we have w. comments in line is below. unfortunately fireSingleImage doesn't have a CRS property so I can't even make that assignment. ''' – William Ross May 21 '20 at 18:07
  • i've included our full code if you'd like to take a look – William Ross May 21 '20 at 18:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.