Lecture 11: Environmental Sciences Packages¶
A collection of Python packages for atmospheric and environmental sciences¶
Statistics¶
- SciPy: SciPy provides algorithms for optimization, integration, interpolation, eigenvalue problems, algebraic equations, differential equations, statistics and many other classes of problems.
- statsmodels: A python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration.
Visualizations¶
- Seaborn: A Python data visualization library based on matplotlib. It provides a high-level interface for drawing attractive and informative statistical graphics.
- Bokeh An interactive visualization library for modern web browsers.
- plotly: An interactive, open-source, and browser-based graphing library for Python.
Machine Learning:¶
- Scikit-learn: A machine learning library that provides almost all the machine learning algorithms you might need.
- PyCaret: PyCaret is an open-source, low-code machine learning library in Python that automates machine learning workflows
- TensorFlow: An end-to-end open source platform for machine learning.TensorFlow is a framework for defining and running computations that involve tensors, which are partially defined computational objects that eventually produce a value.
Geospatial Analysis and Mapping¶
- Geopandas: GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types.
- xESMF: Universal regridder for geospatial data.
- pyResample: Resampling geospatial image data.
- Rasterio: Rasterio is a GDAL and Numpy-based Python library designed to make your work with geospatial raster data more productive, more fun — more Zen. It is a highly useful module for raster processing which you can use for reading and writing several different raster formats in Python. Python automatically registers all known GDAL drivers for reading supported formats when importing the module.
- Regionmask: Plotting and creation of masks of spatial regions
- GeoViews: Library that makes it easy to explore and visualize geographical, meteorological, and oceanographic datasets, such as those used in weather, climate, and remote sensing research
Hydrology:¶
- Pytesmo: Python Toolbox for the Evaluation of Soil Moisture Observations.
- pyDEM: PyDEM is a package for topographic (terrain) analysis. It takes in digital elevation model (DEM) rasters, and it outputs quantities like slope, aspect, upstream area, and topographic wetness index.
- HydrPy: A framework for the development and application of hydrological models based on Python.
- PyFlo: It is an open-source Python library for performing hydraulic and hydrology stormwater analysis. Features include network hydraulic grade analysis and time/iteration based storage and flood routing simulations.
Atmospheric Science:¶
- MetPy: It is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
- ACT: Toolkit for working with atmospheric time-series datasets of varying dimensions.
- Satpy: Reading, manipulating, and writing data from remote-sensing earth-observing meteorological satellite instruments.
- climmetlab: Python package aiming at simplifying access to climate and meteorological datasets
- xgcm: Extends the xarray data model to understand finite volume grid cells (common in General Circulation Models) and provides interpolation and difference operations for such grids
Geology¶
- GemPy: GemPy is a tool for generating three-dimensional structural geological models in Python. It allows the user to create complex combinations of stratigraphical and structural features such as folds, faults, and unconformities. It was furthermore designed to enable probabilistic modeling to address parameter and model uncertainties.
- MTpy: MTpy is a Python Toolbox for Magnetotelluric (MT) Data Processing, Analysis, Modelling and Visualization.
- Pyleoclim: Pyleoclim is a Python package designed for the analysis of paleoclimate data. Pyleoclim leverages various data science libraries (numpy, pandas, scikit-learn) for time series analysis, as well as and Matplotlib and Cartopy for the creation of publication-quality figures.
Geopandas¶
GeoPandas, as the name suggests, extends the popular data science library pandas by adding support for geospatial data.
The core data structure in GeoPandas is the geopandas.GeoDataFrame, a subclass of pandas.DataFrame, that can store geometry columns and perform spatial operations. The geopandas.GeoSeries, a subclass of pandas.Series, handles the geometries. Therefore, your GeoDataFrame is a combination of pandas.Series, with traditional data (numerical, boolean, text etc.), and geopandas.GeoSeries, with geometries (points, polygons etc.). You can have as many columns with geometries as you wish; there's no limit typical for desktop GIS software.
Each GeoSeries can contain any geometry type (you can even mix them within a single array) and has a GeoSeries.crs attribute, which stores information about the projection (CRS stands for Coordinate Reference System). Therefore, each GeoSeries in a GeoDataFrame can be in a different projection, allowing you to have, for example, multiple versions (different projections) of the same geometry.
Only one GeoSeries in a GeoDataFrame is considered the active geometry, which means that all geometric operations applied to a GeoDataFrame operate on this active column.
Reading files¶
In [1]:
Copied!
import geopandas as gpd
gdf = gpd.read_file('/Users/xiaomengjin/Dropbox/0_Rutgers/3_Teaching/Research_Computing/NCA5_Atlas_Global_Warming_Level_1/NCA_Atlas_Counties.shp')
gdf
import geopandas as gpd
gdf = gpd.read_file('/Users/xiaomengjin/Dropbox/0_Rutgers/3_Teaching/Research_Computing/NCA5_Atlas_Global_Warming_Level_1/NCA_Atlas_Counties.shp')
gdf
Out[1]:
| OBJECTID | NAME | STATE_NAME | STATE_ABBR | FIPS | pr_above_n | prmax1day_ | prmax5yr_G | tavg_GWL1 | tmax1day_G | ... | tmean_jja_ | tmin_days_ | tmin_day_1 | tmin_day_2 | tmin_jja_G | pr_annual_ | pr_days_ab | SHAPE_Leng | SHAPE_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | Autauga County | Alabama | AL | 01001 | 8.866368 | 1.204632 | 0.560895 | 1.732316 | 1.543474 | ... | 1.604421 | 22.678737 | -0.039947 | -8.215921 | 1.774947 | 2.231684 | 11.870172 | 2.062534 | 0.150258 | POLYGON ((-86.41312 32.70739, -86.41219 32.526... |
| 1 | 2 | Baldwin County | Alabama | AL | 01003 | 10.447559 | 3.454147 | 1.104216 | 1.633755 | 1.526716 | ... | 1.508863 | 19.807863 | -0.008627 | -5.261471 | 1.667588 | 3.200307 | 11.634376 | 9.150287 | 0.398401 | MULTIPOLYGON (((-87.96018 30.66235, -87.96046 ... |
| 2 | 3 | Barbour County | Alabama | AL | 01005 | 10.930474 | 0.800702 | 0.499965 | 1.649667 | 1.629789 | ... | 1.571404 | 23.813614 | -0.017614 | -6.758772 | 1.625246 | 2.796249 | 16.953608 | 2.681671 | 0.223264 | POLYGON ((-85.25784 32.14794, -85.25924 32.145... |
| 3 | 4 | Bibb County | Alabama | AL | 01007 | 10.956795 | 2.631909 | 1.330318 | 1.789705 | 1.623523 | ... | 1.684205 | 22.352250 | -0.083273 | -9.176136 | 1.856455 | 1.619284 | 14.511721 | 1.887436 | 0.156487 | POLYGON ((-87.06574 33.24691, -87.02685 33.246... |
| 4 | 5 | Blount County | Alabama | AL | 01009 | 18.316114 | 6.689886 | 6.570682 | 1.757045 | 1.530205 | ... | 1.639023 | 18.284318 | -0.132205 | -9.861114 | 1.868750 | 3.069511 | 23.057712 | 2.413198 | 0.164411 | POLYGON ((-86.45302 34.25932, -86.44414 34.259... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3106 | 3107 | Washakie County | Wyoming | WY | 56043 | 17.864115 | 6.655069 | 8.193466 | 2.103356 | 2.002747 | ... | 2.147466 | 0.158741 | -5.406983 | -11.806500 | 2.097448 | 5.085914 | 16.041836 | 4.247988 | 0.650772 | POLYGON ((-108.54727 44.16848, -108.51002 44.1... |
| 3107 | 3108 | Weston County | Wyoming | WY | 56045 | 20.280563 | 5.863631 | 8.519034 | 2.067239 | 2.078767 | ... | 2.164869 | 0.597778 | -4.788625 | -10.939795 | 2.158937 | 4.700654 | 18.228163 | 3.408275 | 0.695663 | POLYGON ((-104.05450 44.18039, -104.05470 43.9... |
| 3108 | 3109 | NaN | Puerto Rico | PR | 72 | 4.800000 | 0.700000 | 5.800000 | 0.950000 | NaN | ... | NaN | 28.700000 | NaN | NaN | NaN | NaN | NaN | 10.736738 | 0.764291 | MULTIPOLYGON (((-66.53130 17.88292, -66.53401 ... |
| 3109 | 3110 | NaN | Alaska | AK | 02 | 22.300000 | 6.500000 | 5.900000 | 2.100000 | NaN | ... | NaN | NaN | -7.449000 | -12.600000 | NaN | NaN | NaN | 557.646894 | 281.718582 | MULTIPOLYGON (((-179.10933 51.30098, -179.1060... |
| 3110 | 3111 | NaN | Hawaii | HI | 15 | 9.200000 | 5.300000 | 8.100000 | 1.070000 | NaN | ... | NaN | 29.600000 | NaN | -1.100000 | NaN | NaN | NaN | 18.946659 | 1.440200 | MULTIPOLYGON (((-155.06851 19.72998, -155.0680... |
3111 rows × 23 columns
In [2]:
Copied!
### Pandas indexing and selection methods also work with geopandas
gdf_nj = gdf.loc[gdf.STATE_NAME == 'New Jersey']
### Pandas indexing and selection methods also work with geopandas
gdf_nj = gdf.loc[gdf.STATE_NAME == 'New Jersey']
In [3]:
Copied!
gdf_nj = gdf_nj.set_index('NAME')
gdf_nj = gdf_nj.set_index('NAME')
In [4]:
Copied!
gdf_nj
gdf_nj
Out[4]:
| OBJECTID | STATE_NAME | STATE_ABBR | FIPS | pr_above_n | prmax1day_ | prmax5yr_G | tavg_GWL1 | tmax1day_G | tmax_days_ | ... | tmean_jja_ | tmin_days_ | tmin_day_1 | tmin_day_2 | tmin_jja_G | pr_annual_ | pr_days_ab | SHAPE_Leng | SHAPE_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NAME | |||||||||||||||||||||
| Atlantic County | 1741 | New Jersey | NJ | 34001 | 12.709769 | 5.173795 | 9.893641 | 1.964590 | 1.760179 | 0.464282 | ... | 1.791410 | 9.996231 | -0.214846 | -14.452641 | 1.862923 | 2.726068 | 15.623388 | 12.040877 | 0.152625 | MULTIPOLYGON (((-74.41362 39.54437, -74.41400 ... |
| Bergen County | 1742 | New Jersey | NJ | 34003 | 16.508688 | 5.163812 | 8.624438 | 2.108313 | 1.919000 | 0.723812 | ... | 1.930625 | 9.349438 | -0.502000 | -14.903813 | 1.896625 | 4.355973 | 13.356645 | 2.186622 | 0.065902 | MULTIPOLYGON (((-74.12046 40.85576, -74.12117 ... |
| Burlington County | 1743 | New Jersey | NJ | 34005 | 17.037017 | 7.749103 | 12.298672 | 2.017741 | 1.757276 | 0.894690 | ... | 1.831379 | 8.895655 | -0.335707 | -14.852569 | 1.884414 | 4.160636 | 20.504922 | 3.030257 | 0.222051 | MULTIPOLYGON (((-74.44213 39.55354, -74.44109 ... |
| Camden County | 1744 | New Jersey | NJ | 34007 | 17.107375 | 7.184875 | 12.570375 | 2.048250 | 1.841625 | 0.987312 | ... | 1.908937 | 10.710687 | -0.242750 | -14.977312 | 1.952188 | 4.326228 | 19.452632 | 1.772704 | 0.061134 | MULTIPOLYGON (((-75.08716 39.97180, -75.08736 ... |
| Cape May County | 1745 | New Jersey | NJ | 34009 | 13.507833 | 5.575222 | 10.194389 | 1.884833 | 1.601944 | 0.269222 | ... | 1.729889 | 10.808944 | -0.200167 | -14.010056 | 1.793889 | 2.755018 | 18.012068 | 9.960371 | 0.069066 | MULTIPOLYGON (((-74.58412 39.30441, -74.58470 ... |
| Cumberland County | 1746 | New Jersey | NJ | 34011 | 16.687061 | 7.617212 | 11.978455 | 1.981485 | 1.765152 | 0.580061 | ... | 1.856576 | 11.125091 | -0.220667 | -14.537727 | 1.934333 | 3.991424 | 21.338093 | 5.241841 | 0.135073 | MULTIPOLYGON (((-75.41831 39.41545, -75.41792 ... |
| Essex County | 1747 | New Jersey | NJ | 34013 | 14.555750 | 5.560125 | 10.055875 | 2.089000 | 1.897375 | 0.731250 | ... | 1.935875 | 8.622875 | -0.493250 | -14.683375 | 1.923375 | 3.768749 | 10.525440 | 1.138257 | 0.035286 | POLYGON ((-74.32281 40.90884, -74.32247 40.908... |
| Gloucester County | 1748 | New Jersey | NJ | 34015 | 16.226667 | 8.225500 | 12.450750 | 2.078042 | 1.969333 | 1.027333 | ... | 2.018917 | 13.232708 | -0.162542 | -15.102375 | 2.071208 | 4.281921 | 20.030672 | 2.053469 | 0.089123 | MULTIPOLYGON (((-75.40662 39.78106, -75.40773 ... |
| Hudson County | 1749 | New Jersey | NJ | 34017 | 15.946667 | 8.062000 | 11.879333 | 2.177333 | 2.002333 | 0.906333 | ... | 2.037333 | 14.763666 | -0.102000 | -15.584333 | 2.029667 | 3.864291 | 13.227471 | 1.524890 | 0.012900 | MULTIPOLYGON (((-74.16091 40.64526, -74.16016 ... |
| Hunterdon County | 1750 | New Jersey | NJ | 34019 | 19.074233 | 7.458533 | 11.135267 | 2.087933 | 1.913300 | 0.621533 | ... | 1.938633 | 4.455567 | -1.075267 | -14.064233 | 1.929167 | 4.933982 | 16.882525 | 1.775757 | 0.120479 | POLYGON ((-74.83809 40.75226, -74.82749 40.744... |
| Mercer County | 1751 | New Jersey | NJ | 34021 | 17.423600 | 7.210800 | 9.612867 | 2.075733 | 1.844733 | 0.839733 | ... | 1.949600 | 8.515200 | -0.354933 | -15.028600 | 1.959933 | 4.484226 | 15.734961 | 1.534344 | 0.062769 | POLYGON ((-74.72324 40.37739, -74.72206 40.375... |
| Middlesex County | 1752 | New Jersey | NJ | 34023 | 16.515458 | 7.180292 | 11.649917 | 2.084208 | 1.860917 | 0.883833 | ... | 1.975458 | 9.029208 | -0.369500 | -14.953375 | 1.981625 | 4.738598 | 13.874158 | 3.026746 | 0.085986 | MULTIPOLYGON (((-74.34076 40.48271, -74.34389 ... |
| Monmouth County | 1753 | New Jersey | NJ | 34025 | 17.070813 | 7.863531 | 13.860406 | 2.005531 | 1.817750 | 0.617469 | ... | 1.890375 | 9.329031 | -0.280687 | -14.637719 | 1.907125 | 3.912146 | 15.619226 | 4.256690 | 0.130309 | MULTIPOLYGON (((-73.98408 40.41741, -73.99008 ... |
| Morris County | 1754 | New Jersey | NJ | 34027 | 17.161424 | 5.485061 | 9.291970 | 2.105303 | 1.951606 | 0.370364 | ... | 1.962485 | 3.653424 | -1.411697 | -13.954576 | 1.922727 | 4.719114 | 14.513126 | 2.380815 | 0.133213 | POLYGON ((-74.50049 41.08601, -74.49997 41.085... |
| Ocean County | 1755 | New Jersey | NJ | 34029 | 15.653977 | 7.068045 | 12.782796 | 1.972886 | 1.729659 | 0.614295 | ... | 1.797545 | 8.254818 | -0.323955 | -14.488704 | 1.836182 | 3.356637 | 17.491542 | 15.648201 | 0.173052 | MULTIPOLYGON (((-74.09657 40.12414, -74.09687 ... |
| Passaic County | 1756 | New Jersey | NJ | 34031 | 15.799063 | 4.963187 | 9.094875 | 2.083125 | 1.851500 | 0.403250 | ... | 1.885688 | 4.598750 | -1.223063 | -13.993937 | 1.824375 | 4.650235 | 13.346869 | 1.556083 | 0.054999 | POLYGON ((-74.23429 41.14302, -74.21777 41.136... |
| Salem County | 1757 | New Jersey | NJ | 34033 | 17.727826 | 9.147435 | 12.954217 | 2.074696 | 1.935739 | 0.937087 | ... | 2.003348 | 13.687652 | -0.174261 | -15.161043 | 2.073043 | 4.786771 | 22.100363 | 4.394849 | 0.092340 | MULTIPOLYGON (((-75.48224 39.65725, -75.48268 ... |
| Somerset County | 1758 | New Jersey | NJ | 34035 | 18.454240 | 6.473880 | 9.478840 | 2.098480 | 1.962760 | 0.856760 | ... | 2.016280 | 6.486240 | -0.787800 | -14.321200 | 2.005080 | 4.908482 | 18.333972 | 1.699680 | 0.083972 | POLYGON ((-74.55469 40.75678, -74.55335 40.756... |
| Sussex County | 1759 | New Jersey | NJ | 34037 | 16.832972 | 4.646583 | 7.895972 | 2.144833 | 1.996778 | 0.282000 | ... | 2.002611 | 2.364194 | -2.251417 | -13.478472 | 1.904444 | 4.451370 | 14.487146 | 1.685282 | 0.148757 | POLYGON ((-74.67081 41.34637, -74.63760 41.331... |
| Union County | 1760 | New Jersey | NJ | 34039 | 14.994167 | 6.106333 | 11.412667 | 2.124167 | 1.841167 | 0.860667 | ... | 1.944667 | 8.997667 | -0.469667 | -15.386333 | 1.940500 | 4.059828 | 12.067153 | 1.369704 | 0.028546 | POLYGON ((-74.36903 40.73927, -74.36647 40.737... |
| Warren County | 1761 | New Jersey | NJ | 34041 | 19.093179 | 6.100571 | 7.801143 | 2.126071 | 2.035357 | 0.390464 | ... | 2.010179 | 2.703071 | -1.738714 | -13.543821 | 1.961036 | 4.861140 | 15.735885 | 1.741006 | 0.100309 | POLYGON ((-74.93889 41.06878, -74.90189 41.034... |
21 rows × 22 columns
Simple accessors and methods¶
Now we have our GeoDataFrame and can start working with its geometry.
Since there was only one geometry column in the NCA dataset, this column automatically becomes the active geometry and spatial methods used on the GeoDataFrame will be applied to the "geometry" column.
In [5]:
Copied!
gdf_nj.geometry
gdf_nj.geometry
Out[5]:
NAME Atlantic County MULTIPOLYGON (((-74.41362 39.54437, -74.41400 ... Bergen County MULTIPOLYGON (((-74.12046 40.85576, -74.12117 ... Burlington County MULTIPOLYGON (((-74.44213 39.55354, -74.44109 ... Camden County MULTIPOLYGON (((-75.08716 39.97180, -75.08736 ... Cape May County MULTIPOLYGON (((-74.58412 39.30441, -74.58470 ... Cumberland County MULTIPOLYGON (((-75.41831 39.41545, -75.41792 ... Essex County POLYGON ((-74.32281 40.90884, -74.32247 40.908... Gloucester County MULTIPOLYGON (((-75.40662 39.78106, -75.40773 ... Hudson County MULTIPOLYGON (((-74.16091 40.64526, -74.16016 ... Hunterdon County POLYGON ((-74.83809 40.75226, -74.82749 40.744... Mercer County POLYGON ((-74.72324 40.37739, -74.72206 40.375... Middlesex County MULTIPOLYGON (((-74.34076 40.48271, -74.34389 ... Monmouth County MULTIPOLYGON (((-73.98408 40.41741, -73.99008 ... Morris County POLYGON ((-74.50049 41.08601, -74.49997 41.085... Ocean County MULTIPOLYGON (((-74.09657 40.12414, -74.09687 ... Passaic County POLYGON ((-74.23429 41.14302, -74.21777 41.136... Salem County MULTIPOLYGON (((-75.48224 39.65725, -75.48268 ... Somerset County POLYGON ((-74.55469 40.75678, -74.55335 40.756... Sussex County POLYGON ((-74.67081 41.34637, -74.63760 41.331... Union County POLYGON ((-74.36903 40.73927, -74.36647 40.737... Warren County POLYGON ((-74.93889 41.06878, -74.90189 41.034... Name: geometry, dtype: geometry
In [6]:
Copied!
gdf_nj.loc['Middlesex County'].geometry
gdf_nj.loc['Middlesex County'].geometry
Out[6]:
Getting polygon boundary and centroid¶
To get the boundary of each polygon (LineString), access the GeoDataFrame.boundary:
In [7]:
Copied!
gdf_nj["boundary"] = gdf_nj.boundary
gdf_nj["boundary"]
gdf_nj["boundary"] = gdf_nj.boundary
gdf_nj["boundary"]
Out[7]:
NAME Atlantic County MULTILINESTRING ((-74.41362 39.54437, -74.4140... Bergen County MULTILINESTRING ((-74.12046 40.85576, -74.1211... Burlington County MULTILINESTRING ((-74.44213 39.55354, -74.4410... Camden County MULTILINESTRING ((-75.08716 39.97180, -75.0873... Cape May County MULTILINESTRING ((-74.58412 39.30441, -74.5847... Cumberland County MULTILINESTRING ((-75.41831 39.41545, -75.4179... Essex County LINESTRING (-74.32281 40.90884, -74.32247 40.9... Gloucester County MULTILINESTRING ((-75.40662 39.78106, -75.4077... Hudson County MULTILINESTRING ((-74.16091 40.64526, -74.1601... Hunterdon County LINESTRING (-74.83809 40.75226, -74.82749 40.7... Mercer County LINESTRING (-74.72324 40.37739, -74.72206 40.3... Middlesex County MULTILINESTRING ((-74.34076 40.48271, -74.3438... Monmouth County MULTILINESTRING ((-73.98408 40.41741, -73.9900... Morris County LINESTRING (-74.50049 41.08601, -74.49997 41.0... Ocean County MULTILINESTRING ((-74.09657 40.12414, -74.0968... Passaic County LINESTRING (-74.23429 41.14302, -74.21777 41.1... Salem County MULTILINESTRING ((-75.48224 39.65725, -75.4826... Somerset County LINESTRING (-74.55469 40.75678, -74.55335 40.7... Sussex County LINESTRING (-74.67081 41.34637, -74.63760 41.3... Union County LINESTRING (-74.36903 40.73927, -74.36647 40.7... Warren County LINESTRING (-74.93889 41.06878, -74.90189 41.0... Name: boundary, dtype: geometry
In [8]:
Copied!
gdf_nj.loc['Middlesex County'].boundary
gdf_nj.loc['Middlesex County'].boundary
Out[8]:
Since we have saved boundary as a new column, we now have two geometry columns in the same GeoDataFrame.
We can also create new geometries, which could be, for example, the centroid:
In [9]:
Copied!
gdf_nj["centroid"] = gdf_nj.centroid
gdf_nj["centroid"]
gdf_nj["centroid"] = gdf_nj.centroid
gdf_nj["centroid"]
/var/folders/7b/6t7qqfj57bb0_ml_y_5bw86r0000gn/T/ipykernel_21619/3860640333.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf_nj["centroid"] = gdf_nj.centroid
Out[9]:
NAME Atlantic County POINT (-74.67735 39.48279) Bergen County POINT (-74.07824 40.96239) Burlington County POINT (-74.66766 39.87828) Camden County POINT (-74.95778 39.80156) Cape May County POINT (-74.80712 39.15175) Cumberland County POINT (-75.11028 39.37396) Essex County POINT (-74.24861 40.78830) Gloucester County POINT (-75.13712 39.71375) Hudson County POINT (-74.07725 40.74073) Hunterdon County POINT (-74.91224 40.56730) Mercer County POINT (-74.70174 40.28344) Middlesex County POINT (-74.41303 40.43831) Monmouth County POINT (-74.22610 40.25884) Morris County POINT (-74.54452 40.86201) Ocean County POINT (-74.29942 39.91042) Passaic County POINT (-74.30083 41.03439) Salem County POINT (-75.34574 39.58701) Somerset County POINT (-74.61633 40.56351) Sussex County POINT (-74.69084 41.13931) Union County POINT (-74.31056 40.66027) Warren County POINT (-74.99735 40.85715) Name: centroid, dtype: geometry
Plotting Data on Maps¶
GeoPandas can also plot maps, so we can check how the geometries appear in space. To plot the active geometry, call GeoDataFrame.plot(). To color code by data, pass in that column as the first argument.
In [10]:
Copied!
gdf_nj.plot()
gdf_nj.plot()
Out[10]:
<Axes: >
In [11]:
Copied!
gdf_nj.plot('tavg_GWL1', legend = True)
gdf_nj.plot('tavg_GWL1', legend = True)
Out[11]:
<Axes: >
In [12]:
Copied!
gdf_nj.columns
gdf_nj.columns
Out[12]:
Index(['OBJECTID', 'STATE_NAME', 'STATE_ABBR', 'FIPS', 'pr_above_n',
'prmax1day_', 'prmax5yr_G', 'tavg_GWL1', 'tmax1day_G', 'tmax_days_',
'tmax_day_1', 'tmax_day_2', 'tmean_jja_', 'tmin_days_', 'tmin_day_1',
'tmin_day_2', 'tmin_jja_G', 'pr_annual_', 'pr_days_ab', 'SHAPE_Leng',
'SHAPE_Area', 'geometry', 'boundary', 'centroid'],
dtype='object')
You can also explore your data interactively using GeoDataFrame.explore(), which behaves in the same way plot() does but returns an interactive map instead.
Switching the active geometry (GeoDataFrame.set_geometry) to centroids, we can plot the same data using point geometry.
In [13]:
Copied!
gdf_nj = gdf_nj.set_geometry('boundary')
gdf_nj = gdf_nj.set_geometry('boundary')
In [14]:
Copied!
gdf_nj.plot('tavg_GWL1', legend = True)
gdf_nj.plot('tavg_GWL1', legend = True)
Out[14]:
<Axes: >
And we can also layer both GeoSeries on top of each other. We just need to use one plot as an axis for the other.
In [15]:
Copied!
ax = gdf_nj["geometry"].plot()
gdf_nj["centroid"].plot(ax=ax, color="black")
ax = gdf_nj["geometry"].plot()
gdf_nj["centroid"].plot(ax=ax, color="black")
Out[15]:
<Axes: >
Now we set the active geometry back to the original GeoSeries.
In [16]:
Copied!
gdf_nj = gdf_nj.set_geometry("geometry")
gdf_nj = gdf_nj.set_geometry("geometry")
Buffer¶
In other cases, we may need to buffer the geometry using GeoDataFrame.buffer(). Geometry methods are automatically applied to the active geometry, but we can apply them directly to any GeoSeries as well. Let's buffer the counties and their centroids and plot both on top of each other.
In [17]:
Copied!
# buffering the active geometry by 0.1 degrees
gdf_nj["buffered"] = gdf_nj.buffer(0.1)
gdf_nj["buffered_centroid"] = gdf_nj["centroid"].buffer(0.05)
# buffering the active geometry by 0.1 degrees
gdf_nj["buffered"] = gdf_nj.buffer(0.1)
gdf_nj["buffered_centroid"] = gdf_nj["centroid"].buffer(0.05)
/var/folders/7b/6t7qqfj57bb0_ml_y_5bw86r0000gn/T/ipykernel_21619/3355730635.py:2: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf_nj["buffered"] = gdf_nj.buffer(0.1) /var/folders/7b/6t7qqfj57bb0_ml_y_5bw86r0000gn/T/ipykernel_21619/3355730635.py:4: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf_nj["buffered_centroid"] = gdf_nj["centroid"].buffer(0.05)
In [18]:
Copied!
# saving the first plot as an axis and setting alpha (transparency) to 0.5
ax = gdf_nj["buffered"].plot(alpha=0.5)
# passing the first plot as an axis to the second
gdf_nj["buffered_centroid"].plot(ax=ax, color="red", alpha=0.5)
# # passing the first plot and setting linewidth to 0.5
gdf_nj["boundary"].plot(ax=ax, color="white", linewidth=0.5)
# saving the first plot as an axis and setting alpha (transparency) to 0.5
ax = gdf_nj["buffered"].plot(alpha=0.5)
# passing the first plot as an axis to the second
gdf_nj["buffered_centroid"].plot(ax=ax, color="red", alpha=0.5)
# # passing the first plot and setting linewidth to 0.5
gdf_nj["boundary"].plot(ax=ax, color="white", linewidth=0.5)
Out[18]:
<Axes: >
Geometry relations¶
We can also ask about the spatial relations of different geometries. Using the geometries above, we can check which of the buffered counties intersect the original geometry of Middlesex, i.e., is within 0.1 degree difference from Middlesex.
First, we get a polygon of Middlesex.
In [19]:
Copied!
middlesex = gdf_nj.loc["Middlesex County", "geometry"]
middlesex
middlesex = gdf_nj.loc["Middlesex County", "geometry"]
middlesex
Out[19]:
The polygon is a shapely geometry object, as any other geometry used in GeoPandas.
In [20]:
Copied!
type(middlesex)
type(middlesex)
Out[20]:
shapely.geometry.multipolygon.MultiPolygon
Then we can check which of the geometries in gdf["buffered"] intersects it.
In [21]:
Copied!
gdf_nj["buffered"].intersects(middlesex)
gdf_nj["buffered"].intersects(middlesex)
Out[21]:
NAME Atlantic County False Bergen County False Burlington County False Camden County False Cape May County False Cumberland County False Essex County True Gloucester County False Hudson County True Hunterdon County False Mercer County True Middlesex County True Monmouth County True Morris County True Ocean County True Passaic County False Salem County False Somerset County True Sussex County False Union County True Warren County False dtype: bool
Integration with Cartopy¶
In [22]:
Copied!
import numpy as np
from matplotlib import pyplot as plt
from cartopy import crs as ccrs
import cartopy
import numpy as np
from matplotlib import pyplot as plt
from cartopy import crs as ccrs
import cartopy
In [24]:
Copied!
central_lat = 37.5
central_lon = -96
extent = [-120, -70, 20, 50.5]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])
crs = ccrs.AlbersEqualArea(central_lon, central_lat)
plt.figure(figsize=(12, 6))
ax = plt.axes(projection=crs)
ax.set_extent(extent)
crs_proj4 = crs.proj4_init
gdf.to_crs(crs_proj4).plot('tavg_GWL1', vmin = 1, vmax = 2.5, ax = ax)
ax.add_feature(cartopy.feature.STATES, edgecolor='black')
ax.coastlines()
central_lat = 37.5
central_lon = -96
extent = [-120, -70, 20, 50.5]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])
crs = ccrs.AlbersEqualArea(central_lon, central_lat)
plt.figure(figsize=(12, 6))
ax = plt.axes(projection=crs)
ax.set_extent(extent)
crs_proj4 = crs.proj4_init
gdf.to_crs(crs_proj4).plot('tavg_GWL1', vmin = 1, vmax = 2.5, ax = ax)
ax.add_feature(cartopy.feature.STATES, edgecolor='black')
ax.coastlines()
Out[24]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x18d3d5880>
In [ ]:
Copied!