Programming for GIS

import geopandas as gpd
import matplotlib.pyplot as plt

# Load shapefile
gdf = gpd.read_file('countries.shp')

# Display first 5 rows
print(gdf.head())
print(f"CRS: {gdf.crs}")
print(f"Features: {len(gdf)}")

# Plot the data
fig, ax = plt.subplots(figsize=(12, 8))
gdf.plot(ax=ax, edgecolor='black', linewidth=0.5)
plt.title('World Countries')
plt.show()

import geopandas as gpd

# Load data
cities = gpd.read_file('cities.shp')
region = gpd.read_file('region.shp')

# Create 50km buffer around cities
cities_buffered = cities.copy()
cities_buffered['geometry'] = cities.buffer(50000)  # in meters

# Clip buffered cities to region
clipped = gpd.clip(cities_buffered, region)

# Spatial join: find which region each city belongs to
joined = gpd.sjoin(cities, region, how='inner', predicate='within')
print(joined[['city_name', 'region_name']])

# Save result
clipped.to_file('cities_buffered_clipped.shp')

import geopandas as gpd

# Load districts shapefile
districts = gpd.read_file('districts.shp')

# Reproject to UTM for accurate area calculation
districts_utm = districts.to_crs(epsg=32643)

# Calculate area in square kilometers
districts_utm['area_sqkm'] = districts_utm.geometry.area / 1e6

# Dissolve districts by state name
states = districts_utm.dissolve(by='state_name', aggfunc='sum')

# Calculate population density
states['pop_density'] = states['population'] / states['area_sqkm']

# Export as GeoJSON
states.to_file('states_analysis.geojson', driver='GeoJSON')

from osgeo import gdal
import numpy as np

# Open raster file
ds = gdal.Open('dem.tif')
band = ds.GetRasterBand(1)

# Get metadata
print(f'Projection: {ds.GetProjection()}')
print(f'Size: {ds.RasterXSize} x {ds.RasterYSize}')
print(f'Bands: {ds.RasterCount}')

# Read data as NumPy array
data = band.ReadAsArray()
print(f'Min elevation: {np.nanmin(data):.1f} m')
print(f'Max elevation: {np.nanmax(data):.1f} m')

# Get geotransform (coordinates)
gt = ds.GetGeoTransform()
print(f'Origin: ({gt[0]}, {gt[3]})')
print(f'Pixel size: {gt[1]} x {abs(gt[5])}')
ds = None

import rasterio
import numpy as np

# Read NIR and Red bands
with rasterio.open('sentinel2_B08_nir.tif') as nir_src:
    nir = nir_src.read(1).astype(float)
    profile = nir_src.profile

with rasterio.open('sentinel2_B04_red.tif') as red_src:
    red = red_src.read(1).astype(float)

# Calculate NDVI
np.seterr(divide='ignore', invalid='ignore')
ndvi = (nir - red) / (nir + red)
ndvi = np.where(np.isfinite(ndvi), ndvi, 0)

# Save NDVI raster
profile.update(dtype=rasterio.float32, count=1)
with rasterio.open('ndvi_output.tif', 'w', **profile) as dst:
    dst.write(ndvi.astype(np.float32), 1)

from shapely.geometry import Point, LineString, Polygon
from shapely.ops import unary_union

# Create geometries
point = Point(77.2090, 28.6139)  # New Delhi
line = LineString([(77.2, 28.6), (72.8, 19.0)])  # Delhi to Mumbai
polygon = Polygon([(76, 28), (78, 28), (78, 26), (76, 26)])

# Operations
buffer = point.buffer(0.5)  # ~50km buffer in degrees
print(f'Contains point: {polygon.contains(point)}')

# Intersection
intersection = buffer.intersection(polygon)
print(f'Intersection area: {intersection.area:.4f}')

# Union of multiple polygons
poly2 = Polygon([(77, 27), (79, 27), (79, 25), (77, 25)])
merged = unary_union([polygon, poly2])

// Define area of interest
var roi = ee.Geometry.Point([77.2090, 28.6139]);

// Load Sentinel-2 collection
var s2 = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterBounds(roi)
  .filterDate('2023-01-01', '2023-12-31')
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));

// Calculate NDVI for each image
var ndviCollection = s2.map(function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
    .rename('NDVI');
  return ndvi.set('system:time_start', 
    image.get('system:time_start'));
});

// Create time series chart
var chart = ui.Chart.image.series(
  ndviCollection, roi, ee.Reducer.mean(), 10)
  .setOptions({title: 'NDVI Time Series 2023'});
print(chart);

// Display on map
var meanNDVI = ndviCollection.mean();
Map.centerObject(roi, 10);
Map.addLayer(meanNDVI, 
  {min: 0, max: 0.8, palette: ['red','yellow','green']},
  'Mean NDVI');

// Load Sentinel-2 composite
var image = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate('2023-06-01', '2023-09-01')
  .filterBounds(roi)
  .median()
  .select(['B2','B3','B4','B5','B6','B7','B8','B11','B12']);

// Training data (FeatureCollection with class labels)
var training = image.sampleRegions({
  collection: trainingPoints,
  properties: ['landcover'],
  scale: 10
});

// Train Random Forest classifier
var classifier = ee.Classifier.smileRandomForest(100)
  .train(training, 'landcover');

// Classify the image
var classified = image.classify(classifier);

// Display results
Map.addLayer(classified, {
  min: 0, max: 5,
  palette: ['green','darkgreen','blue','yellow','red','gray']
}, 'Land Cover Classification');