2020年11月25日星期三

ogrinfo:操作多个shp文件

 https://gdal.org/drivers/vector/shapefile.html

 

 Normally the OGR Shapefile driver treats a whole directory of shapefiles as a dataset, and a single shapefile within that directory as a layer.

 In this case the directory name should be used as the dataset name. 

However, it is also possible to use one of the files (.shp, .shx or .dbf) in a shapefile set as the dataset name, and then it will be treated as a dataset with one layer.

 

 ogrinfo    -sql "ALTER TABLE  ymss_elements_0  ADD    COLUMN id   integer  "   shapefile_dir

2020年11月23日星期一

ogr2ogr :Spatial Joins

 ogr-spatial-join

OGR command line tools accept only 1 input. But we have 2 inputs for the spatial join. An easy way to fix this, is to use a VRT file. A VRT file allows us to specify multiple inputs and pass them to the command-line tool as layers of a single input.

Unzip the input shapefiles in a single folder on your drive. Create a file named input.vrt in the same folder with the following content.

<OGRVRTDataSource>
    <OGRVRTLayer name="boroughs">
        <SrcDataSource>nybb.shp</SrcDataSource>
        <SrcLayer>nybb</SrcLayer>
    </OGRVRTLayer>
    <OGRVRTLayer name="nursinghomes">
        <SrcDataSource>OEM_NursingHomes_001.shp</SrcDataSource>
        <SrcLayer>OEM_NursingHomes_001</SrcLayer>
    </OGRVRTLayer>
</OGRVRTDataSource>

Open the OSGeo4W shell and cd to the directory containing the shapefiles and the vrt file. Run the ogrinfo command to check if the VRT file is correct.

ogrinfo input.vrt

OGR tools can run SQL queries on the input layers. We will use the ST_INTERSECTS function to find all nursing homes that intersect the boundary of a borough and use the SUM function to find the total nursing home capacity of a borough. Run the following command.

ogrinfo -sql "SELECT b.BoroName, sum(n.Capacity) as total_capacity from
boroughs b, nursinghomes n WHERE ST_INTERSECTS(b.geometry, n.geometry) group
by b.BoroName" -dialect SQLITE input.vrt

You can see that in a single command we got the results by doing a spatial join that takes a lot of clicking around in a GIS environment. We can do a reverse spatial join as well. We can join the name of the Borough to each feature of the Nursing Homes layer. Using the ogr2ogr tool we can write out a shapefile from the resulting join. Note that we are adding a geometry column in the SELECT statement which results in a spatial output. Run the following command:

ogr2ogr -sql "SELECT n.Name, n.Capacity, n.geometry, b.BoroName from
boroughs b, nursinghomes n WHERE ST_INTERSECTS(b.geometry, n.geometry)"
-dialect SQLITE output.shp input.vrt

Open the output.shp in a GIS to verify that the new shapefile as attributes joined from the intersecting borough. You can use ogrinfo command to check that as well.

ogrinfo -al output.shp