Created
August 25, 2018 04:38
-
-
Save nishadhka/aeb7a7af831be9ac84c1a08ea8ebb4a1 to your computer and use it in GitHub Desktop.
NDVI shape file creator from Landsat images, using Google earth engine and Python API
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
import fiona | |
from shapely import geometry | |
import ee | |
ee.Initialize() | |
from gee_library import * | |
from rasterio.features import shapes | |
import geopandas as gp | |
import ntpath | |
import glob | |
import json | |
import pandas as pd | |
shapefiles=glob.glob('*.shp') | |
def path_leaf(path): | |
head, tail = ntpath.split(path) | |
return tail or ntpath.basename(head) | |
collection = ee.ImageCollection('LANDSAT/LC8_L1T_TOA_FMASK') | |
wer=[] | |
for shpfiles in shapefiles: | |
with fiona.open(shpfiles, "r") as shpfile: | |
features = [feature["geometry"] for feature in shpfile] | |
filename=(path_leaf(shpfiles)).split('.')[0] | |
aa=geometry.Polygon(features[0]['coordinates'][0]) | |
bb=aa.bounds | |
ll=(bb[0],bb[1]) | |
ur=(bb[2],bb[3]) | |
nps_bounds = bound_geometry(ll,ur) | |
print filename | |
#print(nps_bounds) | |
grid_collection = collection.filterBounds(nps_bounds) | |
tr_cr = grid_collection.filter(ee.Filter.date('2013-01-01', '2016-01-01')) | |
cc_cr = tr_cr.filter(ee.Filter.lt('CLOUD_COVER', 10)) | |
aa=cc_cr.getInfo() | |
we=[] | |
for aas in aa['features']: | |
we.append(aas['properties']['DATE_ACQUIRED']) | |
print "there are "+str(len(we))+" images for "+filename | |
if len(we) >0: | |
print "hello" | |
tiles,img_affine = img_at_region(geCollection=cc_cr,resolution=30,bands=['B4', 'B5','fmask'],geo_bounds=nps_bounds) | |
ndvi=(tiles['B5']-tiles['B4'])/(tiles['B5']+tiles['B4']) | |
a_ndvi = np.where(ndvi<0.5,0,1) | |
fmask=np.where(tiles['fmask']<=1,0,1) | |
b_ndvi=a_ndvi+fmask | |
c_ndvi=b_ndvi.astype(np.uint8) | |
unique, counts = np.unique(c_ndvi, return_counts=True) | |
dattif=dict(zip(unique, counts)) | |
dattif['date']=we[0] | |
dattif['gridname']=filename | |
wer.append(dattif) | |
mask = None | |
resultsU = ({'properties': {'raster_val': v}, 'geometry': s} for i, (s, v) in enumerate(shapes(c_ndvi, mask=mask, transform=img_affine))) | |
geomsU = list(resultsU) | |
dbU = gp.GeoDataFrame.from_features(geomsU) | |
dbU1=dbU[dbU['raster_val']==1] | |
dbU1.to_file(filename+'.shp', driver='ESRI Shapefile') | |
else: | |
print "no suitable image for "+filename | |
grid_collection =[] | |
tr_cr =[] | |
cc_cr =[] | |
#print filename | |
wer1=pd.DataFrame(wer) | |
wer1.to_csv('/home/sunbird/earthengine/ge_stats.csv') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Though I am not a python expert, the logical streams seems great. Wonderful job