texts: QGIS iface.addVectorLayer() results in incorrect field length

texts:
I ambuilding a python plugin for QGIS. As part of the plugin, I do:
iface.addVectorLayer(geoj, self.selectedLayer, 'ogr')
This is basically adding a layer to the project, with
geoj: a geojson string,
self.selectedLayer: the name I want to give it, and
'ogr': the driver.
I get a layer in my project with most field lengths set to zero:

This is a problem, because when I then try to export the loaded layer as a shapefile or something, I get the following error:
*Export to vector file failed.
Error: Feature write errors:
Error converting value ((1:0)) for attribute field photos_captions: String of length 5 exceeds maximum field length (1)
*
And a complete QGIS Crash!
Do I need to specify some more parameters to be able to load the geojson in my project as a valid layer with correctly defined fields?
Here are the relevant sections of code:
 def getGeoJsonFromSelectedLayer(self):

        if self.selectedLayer:
            url = "https://api.fulcrumapp.com/api/v2/query"

            tableSelector = f"SELECT * FROM \"{self.selectedLayer}\""
            querystring = {"q": tableSelector, "format": "geojson", "headers": "false",
                           "metadata": "false", "arrays": "false", "page": "1", "per_page": "20000"}

            headers = {
                "Accept": "application/json",
                "X-ApiToken": self.API_TOKEN
            }

            response = requests.request(
                "GET", url, headers=headers, params=querystring)

            # For Debugging:
            iface.messageBar().pushMessage(response.text)

            self.createLayerFromGeojson(response.text)

        else:
            iface.messageBar().pushMessage("No Layer Selected")

    def createLayerFromGeojson(self, geoj):
        # if there are features in the list
        if len(geoj) > 0:

            # add the layer to the list of layers
            iface.addVectorLayer(geoj, self.selectedLayer, 'ogr')

        else:
            iface.messageBar().pushMessage("No features found in the geoJSON")


texts: Standard Format for Underpass Clearance?

texts:
Do highway departments typically provide underpass clearance data as part of data sets provided to google maps, et. al.?

If so, can anybody point me at the definition of the data elements, file formats, etc.?


texts: How to see classifications of ArcGIS Pro Iso Cluster Unsupervised Classification output raster?

texts: I input a number of raster bands into the Iso Cluster Unsupervised Classification tool and asked for 5 classifications and specified a signature file to be created. The tool ran for a while and then completed. It put a raster in the Table of Contents that was a single solid color. I looked at the signature file and it had 5 classifications. I tried to look for a way to change the symbology of the output raster to see the classifications, but couldn't find it. How can I get a visualization of the classifications in the actual output raster?


texts: What is the purpose of some functions and files in QGIS Python plugins?

texts:
I'm just trying to understand the architecture of a QGIS Plugin created with the QGIS Plugin Builder, as well as the workflow and the linkage among generated files. If possible, I'd prefer a diagram.

It would also be helpful to know what the following functions do in the Python code. I'm new in QGIS Plugin creation.

__init__( self, iface )
tr( self, message )
add_action( self )
initGui( self )
unload( self )
run( self )


texts: Looking for market share of GIS desktop tools

texts:
What are the relative market shares of the various desktop GISes? I'm interested in this within the open-source and non-opensource realm.

One article that I read, is that the market share of ESRI is estimated to be approximately 30% of the global GIS market in 2009 (Daratech, 2011 - dead link unfortunately) but that they have a much higher perceived 'mind-share' (Batty).

I started thinking about this as I wondered which was the most popular open-source GIS program, and assumed it was QGIS but I'm finding it hard to locate any numbers.

To clarify (in response to a comment), I am interested in GIS Desktop tools that allow quantitative analysis (which would exclude Google Earth).


texts: How to export a 3D map into a image file?

texts:
With qgis, it seems to be possible to visualize 3D with two tools : nviz and Qgis2threejs.

But i would like export the result into an image file. Is it possible to do that with qgis?


texts: Get distance from point to 10 nearest points in ArcGIS

texts: I've a layer with the location of different sampled ponds and a layer with all the ponds in the region. I want to calculate the mean distance to the 10 nearest ponds in the region for each of the sampled ponds. I used point distance to get all the distances but there is no option to only get the 10 nearest distances. Because of the amount of data it will be very intensive to manually search for the 10 nearest.


texts: How to load 2 different WFS layers using OpenLayers?

texts:
I've been trying to load 2 WFS layers with the 3857 projection and I'm only able to load one of them, a line vector layer, while the other (point vector) doesn't load at all.

Later I'm going to want to edit both layers in the web app but for now I'm just trying to get them both to load

The code I'm using to load them is:

var geojsonFormat = new ol.format.GeoJSON();

var vectorSource = new ol.source.Vector({
loader: function(extent, resolution, projection) {
  var url1 =    'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=BFTchambers:bft2&' +
    'outputFormat=text/javascript&format_options=callback:loadFeatures' +
    '&bbox=' + extent.join(',');
$.ajax({url: url1, dataType: 'jsonp', jsonp: false})
    .done(function(response) {
        formatWFS = new ol.format.WFS(),
        sourceVector.addFeatures(formatWFS.readFeatures(response))
        });
},
strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
maxZoom: 20
})),
});

var vectorSource2 = new ol.source.Vector({
loader: function(extent, resolution, projection) {
  var url2 =    'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=BFTchambers:chamber&' +
    'outputFormat=text/javascript&format_options=callback:loadFeatures' +
    '&bbox=' + extent.join(',');
$.ajax({url: url2, dataType: 'jsonp', jsonp: false})
    .done(function(response) {
        pointWFS = new ol.format.WFS(),
        sourceVector2.addFeatures(pointWFS.readFeatures(response))
        });
},
strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
maxZoom: 20
})),
});

/**
 * JSONP WFS callback function.
 * @param {Object} response The response object.
*/
window.loadFeatures = function(response) {
console.log("some stuff loaded");
vectorSource.addFeatures(geojsonFormat.readFeatures(response));
vectorSource2.addFeatures(geojsonFormat.readFeatures(response));
};

selectPointerMove = new ol.interaction.Select({
condition: ol.events.condition.pointerMove
});

var layerVector = new ol.layer.Vector({
source: vectorSource,
style: new ol.style.Style({
image: new ol.style.Circle({
stroke: new ol.style.Stroke({
  color: 'rgba(0, 0, 255, 1.0)',
  width: 5
  })
})
})
});

var layerVector2 = new ol.layer.Vector({
source: vectorSource2,
style: new ol.style.Style({
stroke: new ol.style.Stroke({
  color: 'rgba(0, 0, 255, 1.0)',
  width: 2
})
})
});

var raster = new ol.layer.Tile({
source: new ol.source.BingMaps({
imagerySet: 'Road',
key: 'AsBPC3VYKwiu_XhkUWVqfNQGC6upeSXaFwwVyN9SdMEu8znuAtNkymIX9xzRTfs6'
})
});

var map = new ol.Map({
layers: [raster, layerVector, layerVector2],
target: document.getElementById('map'),
view: new ol.View({
center: [-12000,6711790],
zoom: 12
})
});

Only the layerVector loads, not the layerVector2. How to fix the code?

Entire fiddle here https://jsfiddle.net/Luffydude/2c5o284v/


texts: osgeo ogr layer fieds - have become invisible

texts:
We have e python script, using osgeo, that accesses a PostgreSQL database for GIS data. The script have unexpected stopped working, and I have been able to boil it down to:
.GetFieldCount() returns 0, though fields do exist.

I can reproduce the error via this script, which is a modified version of the online-help.

from osgeo import ogr
import sys

databaseServer = "postgrestest1"
port = "5433" # <-- I know default is 5432, bu this server uses 5433
databaseName = "pgv_2017"
databaseUser = "myname"
databasePW = "mypassw"
connString = "PG: host=%s port=%s dbname=%s user=%s password=%s"    %(databaseServer,port,databaseName,databaseUser,databasePW)

lyr_name = "k.jis"

def GetPGLayerFields( lyr_name ):
    conn = ogr.Open(connString)
    print "conn:", conn

    lyr = conn.GetLayer( lyr_name )
    print "lyr:", lyr

    lyrDefn = lyr.GetLayerDefn()
    print "lyrDefn:", lyrDefn

    num_flds = lyrDefn.GetFieldCount()
    print "num_flds:", num_flds

    for i in range(num_flds):
        print lyrDefn.GetFieldDefn( i ).GetName()

    conn.Destroy()

if __name__ == '__main__':

    GetPGLayerFields( lyr_name )

The output looks like this:

C:\Martin\Work>python test_pg_connect.py
conn: <osgeo.ogr.DataSource; proxy of <Swig Object of type 'OGRDataSourceShadow *' at 0x00000000021C0720> >
lyr: <osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x00000000021C0660> >
lyrDefn: <osgeo.ogr.FeatureDefn; proxy of <Swig Object of type 'OGRFeatureDefnShadow *' at 0x0000000004541030> >
num_flds: 0

C:\Martin\Work>

So the program can see the table, but thinks it has no fields?

But as seen below, the same table seen in the db (with the same user credentials) has 2 fields, 'jis_id' and 'geom'.

Can anybody suggest a way to access these fields from python/osgeo?

/Martin


texts:
I am updating an old demo that used the first versions of Cartodb so that is runs with the latest version of cartodb.js.

I am trying to overlay a set of points on the base map and getting the following error:

TypeError: layer should be provided
cartodb.createLayer(b=Object { map_canvas="map_canvas", map="map", user_name="tobagoborn", more...}, d=undefined, e=undefined, f=undefined)cartodb.js (line 23)
initialize()lndngsites1.html (line 56)
onload(event=load )

Here is my code:

<script type="text/javascript">

function initialize()
{
    var map = new L.Map('map_canvas',
    {
        doubleClickZoom:true,zoomControl:true,scrollWheelZoom: false,touchZoom: false,dragging: true, minZoom: 11, maxZoom: 13
    }).setView(new L.LatLng(11.23,  -60.70), 6);

    var tilesUrl = 'http://tobagoborn.com/tiles/tobagobasemap/{z}/{x}/{y}.png';
    L.tileLayer(tilesUrl, {attribution: "Cartodb", maxZoom: 13}).addTo(map);

    //alert(cartodb.VERSION)

    cartodb_leaflet1 = cartodb.createLayer(
    {
        map_canvas: 'map_canvas',
        map: 'map',
        user_name:"tobagoborn",
        type: cartodb,
        sublayers: [{
            table_name: 'lsites',
            query: "SELECT * FROM {{table_name}}",
            interactivity: "landing_sites,number_of_boats,top_species,site_code"
        }]
  }).addTo(map_canvas)

}
</script>


texts: Python geoprocessing tool for ArcGIS online

texts:
I am looking to create a geoprocessing tool using python in arcGIS 10.2 that will be used in the ArcGIS online web app builder.  This will be my first tool for a web application and I would like some information/advice about how to use ArcGIS REST services and hosted feature services in a python geoprocessing tool.  When creating tools for ArcMap I use the

arcpy.GetParameterAsText()

code and then choose a parameter in the script properties window.  Is there a specific line of code to use to get feature services and REST services and is there a specific parameter to use?

Also how do I show the result in a web app?

My goal is to save out the results as a pdf file at a user defined location or as a separate html page.


texts:
I am writing a QGis 2.18 plugin.
I am designing the plugin with QT, using QGis widgets.  I can use QgsFileWidget to open file, but I did not find how to use it to save file.
If I try to define an output filename from my plugin, I read an error saying the file could not be found, which I believe means that the widget is in 'read' mode not in 'save' mode.

How can I set, from my python code, the QgsFileWidget to operate in save mode, and therefore let me define a filename for saving my output?


texts: Why Server Vector Tiles for MapboxGL Client rendered incorrectly

texts:
I am trying to set up a web server in Node.js that serves vector tiles to be displayed in a browser using MapboxGL JS. The data for the vector tiles is stored in a PostGIS database.

I receive a geojson file and process it createing a database with next sentence:

CREATE TABLE IF NOT EXISTS countries (
  table_id SERIAL,
  properties jsonb not null,
  geom geometry(GeometryZ,4326),
  primary key (table_id)
);

My current set up seems to going in the right direction, as I can see vector tiles being loaded and displayed in the browser. However the rendered result is incorrect (this is a screenshot of a section of my map):

What is rendering incorrectly???

Here is the code:

var express = require('express');
var SphericalMercator = require('sphericalmercator');
var mercator = new SphericalMercator({
    size: 256 //tile size
});

const { pool } = require('../postgressql/config');

var app = express();

app.get('/:namelayer/:z/:x/:y.pbf', (req, res, next) => {
    var options = {
        x: parseInt(req.params.x),
        y: parseInt(req.params.y),
        z: parseInt(req.params.z),
        layerName: req.params.namelayer
    };
    const bbox = mercator.bbox(options.x, options.y, options.z, false, '3857');
    const sql = `
        SELECT ST_AsMVT(q, '${options.layerName}', 4096, 'geom')
        FROM (
            SELECT
                table_id, properties,
                ST_AsMVTGeom(
                    geom,
                    ST_MakeEnvelope(${bbox[0]}, ${bbox[1]}, ${bbox[2]}, ${bbox[3]}, 4326),
                    4096,
                    256,
                    false
                ) AS geom
            FROM ${'public.'+options.layerName} c
        ) q;
    `;

    try {
        pool.query(sql, (error, results) => {
            if (error) {
                return res.status(500).json({
                    ok: false,
                    message: error
                });
            }

            const tile = results.rows[0];
            // set the response header content type
            res.setHeader('Content-Type', 'application/x-protobuf');

            // trigger catch if the vector tile has no data, (return a 204)
            if (tile.st_asmvt.length === 0) {
                res.status(204);
            }

            // send the tile!
            res.send(tile.st_asmvt);
        });
    } catch (e) {
        res.status(404).send({
            error: e.toString(),
        });
    }
});

module.exports = app;

I followed using the idea from this sites:

https://medium.com/nyc-planning-digital/using-the-new-mvt-function-in-postgis-75f8addc1d68
https://geoinquiets.github.io/taller-vt/8_postgis_vt/
https://blog.jawg.io/how-to-make-mvt-with-postgis/


texts: How to add Real Attributes to Shapefile using pyqgis?

texts:
I created the Qmarxan plugin using the fTools plugin as an example base. When I create a new shapefile using the QgsVectorFileWriter function I don't have any difficulty specifying a "Real" field with a width of 24 and a precision of 15. However when adding fields to an existing shapefile the addAttribute call to to a vector layer complains if the field is specified with a length greater than 20 and a precision greater that 5.

Can anyone tell me why this limitation exists and how get around it?

Thanks

TSW


texts: How to connect contour lines from two different sources?

texts: How to connect contour lines from two different sources (LiDAR & Photogarametry) in order to obtain continuous and uniform contour lines layer ? how can i interpolate them together ? How can i Snap lines only if they have the same Height value ?


texts: Find point close to points from postgis table

texts:
I have a postgis point table with a large list of cities.

I need to find duplicates in it, but only if they are 10km close each other, which means that is the same city.

For example, I have 5 points called 'Eneby', and suppose that 4 of this are 10Km close each others, and the 5th is 300Km away, what will mean that is another Eneby town. The query in this case should only have to show 4 points.

So, I've made this query, but is not showing the expected result:

Select p1.id, p2.id, * from points as p1, points as p2 where p1.name='Eneby'
and st_dwithin(p1.geom,p2.geom,0.1, true)
and p1.id=p2.id

The result is the same 5 points, no matter what distance (in degrees) I set. I'm trying with different distance values, but is always the same result.

Does anybody know if is there any mistake in this query?

I've been looking through this forum, and found some similar questions, but can't find exactly this problem.


texts: Convert from eastings and northings to latitude and longitude without ArcGIS

texts:
I am not a GIS user, but I was given some shapefiles to work with. I do not have ArcGIS or any other GIS tool, but I do have Matlab with the mapping toolbox and R.

From what I can tell, the data is in eastings/northings in meters relative to some reference point.  According to the metadata, the projection is WGS_1984_Mercator.  There is another parameter that lists GCS_WGS_1984.  I need to convert this into latitude and longitude relative to zero latitude and longitude.  I do not need very high accuracy - within a degree or two is fine.

Here's some additional info, though I'm not totally sure that it's all relevant:

dbfdata: {[1]  [7.5947e-05]  [0.9000]}
dfbdate: [2012 7 2]
SPHEROID['WGS_1984',6378137.0,298.257223563]]
PROJECTION['Mercator'],
PARAMETER['False_Easting',0.0],
PARAMETER['False_Northing',0.0]

My apologies if I haven't provided enough information - please let me know what else I need to provide.


texts: Wrong Minimum Pixel value using GetStatistics function

texts:
I want to print the minimum and maximum pixel values of a raster.

When I load the raster Pioneertown.row1.asc in Layers Panel in QGIS, it is showing minimum value = 'nan' and maximum value = '1' .

ds = gdal.Open("E:/Sreeraj/Task Global Data CONVERSION/Split/output/Pioneertown_Agriculture/Pioneertown_row1.asc")
stats =   ds.GetRasterBand(1).GetStatistics(0,1)
print stats
[1.0, 1.0, 1.0, 0.0]

According to this Python program in Python Console in QGIS, I am getting minimum value = '1.0', maximum value = '1.0', mean value = '1.0' and standard deviation = '0.0'.

But, I want to get minimum value = 'nan' or '0.0' or '-999' (no data value).

For that, what should I do ?


texts: Defining layer style for KML layer within OGR Virtual Format?

texts:
I have the following OGR virtual format data source, which includes a few GeoJSON files and one KML file:

<OGRVRTDataSource>
  <OGRVRTUnionLayer name="union">

    <GeometryType>wkbGeometryCollection</GeometryType>
    <LayerSRS>EPSG:4326</LayerSRS>

    <OGRVRTLayer name="AMDT_02_2014_Finland_FIR">
      <SrcDataSource>AMDT_02_2014_Finland_FIR.kml</SrcDataSource>
      <SrcSQL>SELECT *,'PEN(c:#FF00FF,w:10px)' AS OGR_STYLE FROM AMDT_02_2014_Finland_FIR</SrcSQL>
    </OGRVRTLayer>

    <OGRVRTLayer name="OGRGeoJSON">
      <SrcDataSource>tsa2014.js</SrcDataSource>
      <SrcSQL>SELECT *,'PEN(c:#00FFFF,w:10px)' AS OGR_STYLE FROM OGRGeoJSON</SrcSQL>
    </OGRVRTLayer>

    <OGRVRTLayer name="OGRGeoJSON">
      <SrcDataSource>ctr2014.js</SrcDataSource>
      <SrcSQL>SELECT *,'PEN(c:#FFFF00,w:5px)' AS OGR_STYLE FROM OGRGeoJSON</SrcSQL>
    </OGRVRTLayer>

  </OGRVRTUnionLayer>
</OGRVRTDataSource>

Setting the layer style for the GeoJSON layers works fine when doing a ogr2ogr conversion to another format such as PDF. However, trying to set the KML layer style with <SrcSQL>SELECT *,'PEN(c:#FF00FF,w:10px)' AS OGR_STYLE FROM AMDT_02_2014_Finland_FIR</SrcSQL> seems to have no effect. Default black pen is used to draw the polygon outline:

Is there a way to change the style of a KML layer within OGR VRT file? At least the syntax above does not seem to work.

The KML file used in this question can be downloaded from Dropbox.


texts: Getting ArcPy autocomplete to work in IDLE?

texts:
I want to write scripts in IDLE, but I love how the built in python window in ArcGIS 10.2 automatically fills in the module names as you start typing them.

For example, if I type arcpy.Lis..., the python window in ArcGIS will suggest arcpy.ListDataStoreItems.

In IDLE, that doesn't occur, and I feel that if I type one wrong letter or mess up a case, it will cause the script to crash and burn.

At the top of my script, I write:

import arcpy.mapping as mapping
# etc

However I don't believe this runs until the script runs, so how can the window know what I've imported. Can I do something in the shell so that the autocomplete works in IDLE?


texts: Position Point Along line at varying percentage?

texts:
Using code I found on Creating point at midpoint of polyline with ArcPy?, I've gotten my code to generate midpoints due to the 0.50 parameter in the position along line tool (below). BUT. I don't want to generate a specific percentage every time. I want to call the percentage in "duration_fraction" associated with each polyline (below), and put it in as the percentage in the points. So Polyline 1 would generate a point 86% of the way along the line. How should I go about this? Can I use a search cursor within a search cursor?

in_data = "L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\Detection1_lines_UTM"
detectpoint = "L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\MEOWWW"
polyline = "L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\Detection1_lines_UTM"

#housekeeping
arcpy.DeleteFeatures_management(detectpoint)
arcpy.env.overwriteOutput = True

#generating the mid point
with arcpy.da.SearchCursor(polyline, "SHAPE@") as in_cursor, \
arcpy.da.InsertCursor(detectpoint, "SHAPE@") as out_cursor:
    for row in in_cursor:
        fraction = row.getValue("duration_fraction")
        midpoint = row[0].positionAlongLine(0.50,True).firstPoint
        out_cursor.insertRow([midpoint])


texts: Rotating the coordinates of a layer by 90 degrees

texts:
I have got a dwg file with the plans for a building, with rooms and everything. (There's been some samples taken inside the building which I want to display.) Now, since it's of a building, the person making this dwg file didn't care much for coordinate systems, and it's all in some (supposedly) arbitrary system centered around the building. This is not an issue.

The problem is that this file (which I have converted to shapefiles) does not honor North being upwards, and that is what I now need to fix. Basically I want all my shapefiles to be rotated 90 degrees, so that North is upwards. It doesn't really help to rotate my data frame, since that would only change North inside my document and not for my shapefiles' coordinates.

Can I change the content of the prj-file in some way, or how could I achieve this?

Using ArcMap 10.2 basic, Python 2.7, FME 2014, Win 7.


texts: Polygon and raster in the same plot do not align

texts:
I can't manage to plot a raster and boundary of country (sp data frame) in the same plot. The data (lon,lat,value) are in mat data frame. The grid is not regular and comes from a climate change ncdf file.

Here is the ncdump original file :

7 variables (excluding dimension variables):
        double lat[rlon,rlat]
            bounds: lat_vertices
            long_name: latitude coordinate
            standard_name: latitude
            units: degrees_north
        double lat_vertices[vertices,rlon,rlat]
            units: degrees_north
        double lon[rlon,rlat]
            bounds: lon_vertices
            long_name: longitude coordinate
            standard_name: longitude
            units: degrees_east
        double lon_vertices[vertices,rlon,rlat]
            units: degrees_east
        double time_bnds[bnds,time]
        float hfss[rlon,rlat,time]
            standard_name: surface_upward_sensible_heat_flux
            long_name: Surface Upward Sensible Heat Flux
            units: W m-2
            original_name: HFX
            cell_methods: time: mean
            missing_value: 1.00000002004088e+20
            _FillValue: 1.00000002004088e+20

            grid_mapping: rotated_latitude_longitude
            coordinates: lat lon
            positive: up
        int rotated_pole[]
            grid_mapping_name: rotated_latitude_longitude
            grid_north_pole_latitude: 39.25
            grid_north_pole_longitude: -162

I try to plot this grid with raster and levelplot package.

library(raster)
library(rasterVis)
library(rworldmap)

data(wrld_simpl)
colnames(mat) <-c("lon","lat","value")

x.min=min(mat$lon)
x.max=max(mat$lon)
y.min=min(mat$lat)
y.max=max(mat$lat)

v.min=min(mat$value)
v.max=max(mat$value)

breaks=seq(v.min,v.max,2)
colour=rev(rainbow(length(breaks), alpha=1))

ras <- raster(nrows=453, ncols=453,xmn=x.min,xmx=x.max, ymn=y.min, ymx=y.max)
rvm <- matrix(mat[,3],453,453,byrow=TRUE)
ras[]=rvm[nrow(rvm):1,]

levelplot

levelplot(ras,at=breaks,col.regions=colour, margin = FALSE) +
layer(sp.polygons(wrld_simpl, fill = NA, col = "grey50"))

and plot

Of course, SP polygon (in WGS84) does not fit to the original ncdf grid. The original data frame is available here :

[https://drive.google.com/file/d/0Byd6YJGGADrfN255ZGhhTjRSUkU/view?usp=sharing][2]


texts: Use of internal weights for Weighted Overlay suitability modeling

texts:
I have derived 'external criterian weights/overall factor weights' applying Pairwise Comparison table. Similarly, pairwise comparison table was used to calculate 'internal weights' separately for each factor.
E.g. for a given suitability modeling, overall factor weight for Landcover resulted as 34% and internal weights were calculated as;
Barren = 41%
Forest = 7%
Low Veg = 39%
Built = 8%
Agriculture land =5%
Now, as listed above internal weights are calculated in percent unit, I am unsure how would I convert them to make compatible with given 1-9 evaluation scale in 'Weighted Overlay' tool enabling me to assign internal weights to each class?


texts:
I have a working QGIS-Server in one of my Docker container. If i use a display with a high resolution (beyond 2000px in width), QGIS-Server is splitting the request into two Requests. The first request with a WIDTH of 2000 and a second request with the rest. This behavior sabotage my use case. I found the setting <WMSMaxWidth type="int"> </WMSMaxWidth>, which i have set to 3000. Unfortunately the result is the same.

Can someone give me a hint or a concrete solution to increase the value in WIDTH, so that everything can be processed with only one WMS GetMap request.


texts: Use field calculator to check for intersection

texts:
Ultimate goal
I plan to retrieve attributes from layers using the field calculator where they intersect with a geometry.
Current step
Checking whether the current $geometry intersects with the layer of interest.
What I tried
intersects($geometry,collect_geometries('Layer')) - and variations of map_get. Since I just learned that intersects has to be handed known geometries, while I want to check IF there is a geometry on the layer, I actually am at a loss. I simply don't know how to check for intersections using the field calculator.
If this is not possible using the field calculator, then maybe Python?


texts:
I would like to know if there is any way of mirroring layers on QGIS. My current project is a bit hard to explain... I'm excavating the inside of a large pot and, in the end I will have a set of plans that correspond to the successive layers I'm exposing as I go along with the excavation. However, for pratical reasons, I Have to do this work upside-down (i.e., the excavation started from the bottom of the pot).

All my records will later be edited in QGIS so, in order to make more understandable for other people, I will need to reverse them. Not just reversing their order, but also flip them horizontally, getting a mirrored image of each layer's plan.

I know I could do it with image processing software, but I was wondering if there is any way of doing it within the QGIS environment - this would probably save me image processing time and make it a more streamlined workflow.

Does anyone has any ideas on this?


texts:
I am trying to implement single_source_dijkstra_path algorithm of the networkx to define the shortest path from a node to all nodes in the network. There are six sources in total, and I want to measure which one of the six nodes has the shortest path to a target node (e.g from node 1 - 6, which node has the shortest distance to node 30). While it works, I think some of the result didn't show the shortest path.

for example, the picture above shows the shortest path to AS Calau, which is to be from Stralsund with 353.3 km. However, with only a gap of a node with a distance of 10 km, the shortest path changes as follow:

From Calau to Bronkow is only 10 km and my network is proven to be all connected (len(list(nx.connected_components)) = 1 ). In Figure 2, it can be seen that the shortest path is no longer from Stralsund but from Rostock Oberseehafen, and the shortest distance is then 641 km. Further analysis also shown that from Stralsund to Bronkow needs to undergo similar path with figure 2, thus resulting in a longer distance.

Does anybody know what I might be doing wrong in this case or does anyone have any explanation for this?


texts: Has anyone tried implementing Project Repair for customised ArcGIS Geoprocessing Tools?

texts:
I'm currently testing a model with a large number of elements (about 400) which takes five input feature classes and three workspaces.

When I test it against my test data locations I can open up the tool dialog and it finds all my data in the test locations that I have set as defaults so quite quickly I can just hit OK to get the tool running.

However, when I move my test data to a new location to simulate what a user of the tool is likely to need to do, then it takes about 30 seconds to locate the data for each of the eight input fields i.e. 4 minutes in total with me having to wait in between.

What I would like to do is to locate the first dataset, and have my tool do something akin to project repair in an MXD which seems to assume that if one dataset moved from disk location A to disk location B then probably the rest did too.

Before I start to investigate writing Python code for tool validation that attempts to do this I am wondering if anybody else has already done so?


texts: Raster DEM vertical units are defined in meters, yet Arcmap displays them in feet

texts:
Background:

I have received a single-band, 32-bit floating-point tif raster DEM, derived from a lidar project.

The tif's XY coordinate system linear units are feet.  Its vertical coordinate system units are in meters.

Here's an Arc desktop screenshot showing those units:

My Question:

When I add this raster to a freshly-opened Arcmap session (and thus no previous CRS), clicking on the raster image with the Identify button reports elevation values in feet.  I was expecting meters.

Additionally, Layer Properties > Statistics reports values in feet, as shown below (I am very familiar with this area, and although the units are not stated, these elevation values are in feet):

A review of the Arc Help documentation sheds no light on this.

Why does Arcmap report elevation values in feet when the vertical units are defined in meters?


texts: Exporting Spatial Data Frame of GWR results from R

texts:
To export Spatial Data Frame of GWR results from R do I need to use readShapeSpatial?

bg<-readShapeSpatial("BG_Div_NO2_Solar_Nozero.shp")
df.bg <- as.data.frame(bg)
df.bg.spdf <- SpatialPolygonsDataFrame(bg, df.bg)

gwr.res1 <- gwr.basic(bg$DIVINDX10 ~ bg$SOLAR_AREA + bg$Avg_NO2_co,
                      data=df.bg.spdf,
                      bw=16,
                      kernel='gaussian',
                      adaptive=TRUE
                      )

readShapeSpatial...


texts: What's the difference between coordinateDimension and spatialDimension?

texts:
The OpenGIS standard mentions three methods to query the dimensions of a Geometry:

dimension(): the inherent dimension of the geometry:

0 for Point
1 for Curve
2 for Surface

coordinateDimension(): the number of coordinates of the geometry:

2 for X,Y
3 for X,Y,Z or X,Y,M
4 for X,Y,Z,M

spatialDimension(): this one has no description in this document.

This RDF document on the opengis.net website further describes coordinateDimension as:

  The number of measurements or axes needed to describe the position of this geometry in a coordinate system.

And describes spatialDimension as:

  The number of measurements or axes needed to describe the spatial position of this geometry in a coordinate system.

So the latter refers to the "spatial position" as opposed to the "position", which doesn't help me much to understand the difference between them.

What is the difference between coordinateDimension() and spatialDimension()?


texts:
Using this implementation of parallelization raises an ImportError: cannot import name 'Parallel' from 'multiprocessing'  The following code tries parallelization with the "denominator" function and should give me the sum of the fields "basalareap","basalareas","basalaread" in a new column. When I import the whole library viafrom multiprocessing import * The process start but comes to no end.

What is wrong with my syntax?

import numpy as np
from multiprocessing import cpu_count, Parallel
import pandas as pd

#Some example dataframe
np.random.seed(4)
layer = pd.DataFrame(np.random.randint(0,25,size=(10, 4)),
                  columns=list(['basalareap','notofinterest', 'basalareas', 'basalaread']))

###Filter Fields by selecting columns of interest
fields = ["basalareap","basalareas","basalaread"]

#In reality data is a geodatframe, it would be:
#layer = layer[fields+['geometry']]

#but here:
layer = fields

data = layer

def denom():
    data['denominator'] = data[["basalareap","basalareas","basalaread"]].sum(axis=1)

cores = cpu_count()
partitions = cores

def parallelize(data,func):
    data_split = np.array_split(data,partitions)
    pool = Pool(cores)
    data = pd.concat(pool.map(func, data_split))
    pool.close()
    pool.join()
    return data

data = parallelize(data,denom)

I'm using Windows 10 and Python 3.7.4


texts: Dissolve multiple shapefiles so that touching polygons become one: GDAL/OGR

texts:
Trying to polygonize large rasters* but since the process is taking forever/crashing in GDAL, I am resorting to doing so in ArcGIS. However there is a 2GB limit on shapefiles in Arc so I am taking a combination approach where I split the rasters into 4 parts with gdalwarp, polygonize them in ArcGIS (polygonizing the quarters in GDAL is also taking forever while the process was extremely fast in Arc), and then putting them back together in GDAL to get around the slowness of gdal_polygonize.py and also the 2GB ArcGIS limit.

I am basing my solution partially on the answer found here, but am running into a problem with the final Dissolve step.

*The pixels of the rasters I am polygonzing only contain 1's (i.e. the pixel value of the things I'm interested in as polygons) and NoData (do not wish to have polygons in these areas), so the attribute table after polygonization looks like this:

Which is want I expect.

The next step is to merge these (4) polygons back together in GDAL like this:

ogr2ogr -f "ESRI Shapefile" merge.shp seg7_9__h06v01_1_2011to2015_min_a.arcpoly.shp -nln merge # create merge.shp
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_b.arcpoly.shp -nln merge # add b.shp to merge.shp
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_c.arcpoly.shp -nln merge
ogr2ogr -f "ESRI Shapefile" -update -append merge.shp seg7_9__h06v01_1_2011to2015_min_d.arcpoly.shp -nln merge

Which works fine and gives me a result that looks like this, where the highlighted polygon is from one of the input shapefile parts and shows the edge between the two:

My last step/goal is to dissolve the merge.shp so that the polygons on the edge that are touching (such as the one highlighted in red and the portion above it) become one polygon. In the post linked above, they suggested using the following command (command1):

ogr2ogr "output.shp" "input.shp" -dialect sqlite -sql "SELECT ST_Union(geometry), field FROM input GROUP BY field"

Because I do not really understand SQL or SQLITE, I wanted to try this command (command2) at first just to see what would happen:

ogr2ogr out_dissolve1.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry), field FROM merge"

but got the following error:

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry), field FROM merge):
 wrong number of arguments to function ST_Union()

I looked online and tried looking at documentation of ST_Union but could not find anything useful.

In conclusion:

Can anyone explain what the SQL statement of the function from command1 referenced above does and
Does anyone know how I can adjust the ogr2ogr command to obtain the desired result, which is dissolving the touching/intersecting polygons into one?

UPDATE: A couple of other variations I've tried and the errors:

I thought using the GRIDCODE field would do the same thing I want because every polygon has a GRIDCODE=1

ogr2ogr dissolve.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry), GRIDCODE FROM merge GROUP BY GRIDCODE"

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry), GRIDCODE FROM merge GROUP BY GRIDCODE):
  wrong number of arguments to function ST_Union()

and

ogr2ogr dissolve.shp merge.shp -dialect sqlite -sql "SELECT ST_Union(geometry) FROM merge"

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry) FROM merge):
  wrong number of arguments to function ST_Union()

Not sure why I'm getting this error or how to fix it.

UPDATE2: I think I found the right command to dissolve based solely on geometry, but I'm getting an error with the ST_Union command:

ogr2ogr dissolve.shp merge.shp -dialect sqlite3 -sql "SELECT ST_Union(geometry) as geometry FROM merge"

ERROR 1: Undefined function 'ST_Union' used.

I tried using GUnion and got the same error.


texts: Putting element with same symbol for changing value in legend with ArcPy?

texts:
I have a layer with forest districts and I need to generate map of every district. That is quite uncomplicated, however maps need to have a legend, where the symbology is the same for every district, but the name changes. Few solutions came to my mind, but they have certain disadvantages:

creating layers for every district and updating symbology, but that seems too be overly complicating task
I could put a text element on the legend and just update its content with name, but that is "not elegant" - eventual moving legend will mess up the layout
I thought, as I use Definition Query anyway, I could use symbology type Single Symbol and just update description with district name, but I couldn`t find entry for this feature.

Here is a sample for second idea:

with arcpy.da.SearchCursor("districts", ["district_name", "SHAPE@"]) as cursor:
    for row in cursor:
        for TextElement in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
            if TextElement.name == "forest_district":
                TextElement.text = row[0]

It would be best if I could treat single symbol description in a similar manner as TextElement, but as I said I couldn't find any entry for this.

So, I would like to know if anyone knows how to edit description in single symbol symbology type with Python? Or maybe someone has better idea for solving my task?

EDIT: I finally put TextElement next to symbol on legend and used the loop from example to assign district names. The good thing is that text can be changed even after grouping legend with TextElement, so there is no mess after moving legend. It looks good, nevertheless I am not very proud of this solution.


texts: How to import images bucket from AWS S3 using python

texts:
How to import images bucket from AWS s3 using python without declaring LOCAL_PATH, i want total bucket images in python environment in array format, from one week on words i am facing this issue, please any one can help out of this.

Here is my python code

import cv2
import numpy as np
import os
from random import shuffle
from tqdm import tqdm
import tensorflow as tf
import matplotlib.pyplot as plt
import tflearn
from tflearn.layers.conv import conv_2d, max_pool_2d
from tflearn.layers.core import input_data, dropout, fully_connected
from tflearn.layers.estimator import regression
import boto
import boto.s3.connection
import sys
import traceback
conn =   boto.connect_s3('XXXXXXXXXXXXXXXX','XXXXXXXXXXXXXXXXXXXXXXXXX')
bucket = conn.get_bucket('dp-xray')

LOCAL_PATH='E:/'

for l in bucket_train:
    keyString = str(l.key)
    d = LOCAL_PATH + keyString
    try:
        l.get_contents_to_filename(d)
    except OSError:
        if not os.path.exists(d):
            os.mkdir(d)

def create_label(image_name):
    """ Create an one-hot encoded vector from image name """
    word_label = image_name.split('.')[-3]
    if word_label == 'Nor':
        return np.array([1,0])
elif word_label == 'Nod':
        return np.array([0,1])

def create_train_data():
training_data = []
for img in tqdm(os.listdir(TRAIN_DIR)):
    path = os.path.join(TRAIN_DIR, img)
    img_data = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
    img_data = cv2.resize(img_data, (IMG_SIZE, IMG_SIZE))
    training_data.append([np.array(img_data), create_label(img)])
shuffle(training_data)
np.save('train_data.npy', training_data)
return training_data

train_data = create_train_data()


texts: How to convert int type OBJECTID to ObjectID type recognised by ArcSDE?

texts: I have an ordinary table with OBJECTID column (type of int) in a MS SQL Server 2008 R2 instance. I need to store it in an ArcSDE geodatabase instance running on top of the SQL Server. The thing is how to make the OBJECTID field to the esri-compliant Object ID type so that feature operations through Feature Server can be applied down the track?


texts: Errors Clipping a Vector: Invalid Geometry

texts:
I have a vector that was polygonized from a raster. I'm hoping to clip it within the boundaries of a certain shapefile "noFertZone" (the turquoise box). In the end I want a "noPKS" shapefile that only encompasses everything within the turquoise box.

However, when I use the clip by mask tool It only clips the elements fully within the turquoise box.

enter image description here

When I do the simple clip tool, like I would in ArcGIS, my input looks like this. However, the clip tool is spitting back an 'invalid geometry' error.

Any ideas on what I'm doing wrong? I've tried every clipping option I could think of.


texts: I have a DEM surface in WGS 1984.I want to calculate surface area and volume in different elevations.I wonder for example If I want to measure the area and volume below some height how shall I set the Plane Height and Z factor to get the results in square meters and cubic meters?


texts:
Is there in Rasterio any way to create a mosaic?

I tried to google it, but I couldn't find anything, maybe someone else already used it.

I want to use some open source possibilities, I already know of possibility to use gdalmerge, gdalwarp.


texts:
This is probably one of the simplest GIS tasks to do, but I am struggling. I have tried to adapt this answer but it appears not to work Create buffers around points and merge with SpatialPolygonsDataFrame to create a list of matched polygons and the area of intersection

Here is my data with long and lat in WGS84

S10.df<-structure(list(id = c(273, 1028, 3635, 6693, 13430), long_wgs = c(-1.564636192,
-1.56677146, -1.554945908, -1.56773222, -1.564188582), lat_wgs = c(53.82838137,
53.82924493, 53.83130225, 53.82848279, 53.82857442)), class = "data.frame", row.names = c("273",
"1023", "3357", "6326", "12301"))

and this code:

stores.sf <- st_as_sf(S10.df, coords = c("long_wgs", "lat_wgs"), crs = 27700, dim = "XY")
stores.buffer <- st_buffer(stores.sf, 500)

to create 500m polygon buffers around each location. When I try to plot I just get one large circle. Further inspection shows:

> stores.buffer$geometry
Geometry set for 5 features
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -501.5677 ymin: -446.1716 xmax: 498.4451 ymax: 553.8313
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
POLYGON ((498.4354 53.82838, 497.7501 27.6604, ...
POLYGON ((498.4332 53.82924, 497.748 27.66127, ...
POLYGON ((498.4451 53.8313, 497.7598 27.66332, ...
POLYGON ((498.4323 53.82848, 497.747 27.6605, 4...
POLYGON ((498.4358 53.82857, 497.7506 27.6606, ...

Which looks an odd set of POLYGON co-ordinates. Has the "transformation" worked so that st_buffer knows the CRS 27700 is in units of meters? Ideally I'd like the "answer" back in WGS84.


texts:
I am using [QGIS ver 2.14.8], Windows 7

I am doing repeated work to check the point feature's attribute .

Clicking one feature and Ctrl + j ( zoom map to the selected rows)

and checking environment from other satellite map layers (to make sure filled information is correct manually)

The points are far away from to each other and the satellite map takes much time to rendering, so I should click or press Ctrl + j to saving time after checking a point .

I can change the shortcut in map window, but it's not possible in attribute table window shortcut


texts: Leaflet inbuild AJAX methods for WMTS getCapabilities request

texts: I was wondering whether the Leaflet library has inbuilt methods to make AJAX requests. My goal is to create a plugin that makes an XMLHttpRequest request to a specified WMTS resource and parse the XML response to JSON notation. If the library does not integrate AJAX methods I would use vanilla JS to require the the information.


texts: Python GDAL package missing header file when installing via pip

texts:
I asked this question on Stack Overflow before I realised this site existed. As this site it more appropriate, I'll ask it here too.

I'm trying to install gdal from pip pip install gdal inside a virtual environment (Ubuntu). It fails because it cannot find cpl_port.h

extensions/gdal_wrap.cpp:2853:22: fatal error: cpl_port.h: No such file or directory
compilation terminated

However GDAL is installed correctly and the header file is located at /usr/include/gdal/cpl_port.h. Is there some environment variable for GDAL that needs to be set in order for pip to find the header files?


texts: Crop multiple rasters using map or mapply?

texts:
I have a folder of rasters with slightly different extents. I would like to apply crop to the whole list of rasters. I'm trying to teach myself how to use the purrr map family, but I can't get that or mapply to apply a function with multiple arguments to a list of files, i.e.

# reproducible example raster w/ categorical values
# from: https://stackoverflow.com/questions/9940495/plot-raster-factir-values-with-ggplot
library(raster)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r <- reclassify(r, c(0, 500, 1,  # from, to, becomes
                     500, 2000, 2))
levels(r)
rcat <- as.factor(r)
levels(rcat)  # 2 levels
extent(rcat)
# make second raster of slightly different extent
extnew1 <- as(extent(178000, 181520, 329400, 334000), "SpatialPolygons")
rcat2 <- crop(rcat, extnew1)
extent(rcat2)

# list of rasters with different extents
mylist <- list(rcat, rcat2)
unique(map(mylist, extent))  # different extents

# crop both to smaller extent
extnew2 <- as(extent(178400, 181500, 329300, 334000), "SpatialPolygons")

# does not work
rcat_crop <- map(mylist, crop, extnew2)

unique(map(rcat_crop, extent)) # matches smallest extent of rasters, not new extent

# does not work
rcat_crop <- mapply(FUN=crop, MoreArgs=list(x=mylist, y=extnew2))

length(rcat_crop)  # =0

What am I doing wrong?

Related question on how to crop using R


texts: How do I create blocks from a road network?

texts: I have a road network and for every enclosed area within the network I want to create a polygon that represents that area. How does one do this?


texts: KeyError: u'M' in Average Nearest Neighbour test

texts:
When I try and run the Average Nearest Neighbour test on my Shapefile I get the error:

KeyError: u'M'
Failed to Execute (AverageNearestNeighbour)

On a "normal" Shapefile this works fine. However I'm trying to use a Shapefile created originally from a CSV with X, Y coordinates.  So the overall workflow I have followed is:

Check format, spacing etc of CSV file. No spaces. Coordinates provided as Eastings and Northings. Made Eastings and Northings into Number format.
In ArcGIS 10.1 - File>Add Data>Add XY Data.
Add the CSV file using Easting and Northing as the X and Y fields. Coordinate system is the same. Click okay.
Right click on the file created >Data>Export Data and save as Shapefile.
Open the Shapefile in ArcGIS, go to Spatial Statistics>Analyzing Patterns>Average Nearest Neighbour, select the new shapefile, click generate report, click okay...and then the error message comes up.

I can't figure out what I'm doing wrong, does anyone have any ideas?

I'm struggling to find out what that error actually means. I've tried deleting all the unnecessary fields from the CSV file (dates etc) in case that was confusing things, but no luck.

I am a relative beginner.

There are 13,600 fields, which I thought might be the problem. But I tried with a vastly reduced number of rows - down to 800 - but it made no difference.


texts: How to do warping of shapefile in GeoTools

texts:
Can any body give me simple warping example in GeoTools?

I have non georeferenced shapefile with polygons and I need to georeference this shapefile by using warp technique.

Please help me.


texts:
I am not sure if this is the right place to ask, if so please recommend a stackexchange,

I am having a lot of issues getting the data approximations that I need for a report.

I am looking for the number of american households where the householder is aged between age 18-35 with annual income between 20,000 and 50,000. And also those with income more than 50,000.


texts: r.univar - which format for working in GRASS GIS?

texts:
As part of a cut and fill-procedure,I want to run the command r.univar, but it returns with no output. I´ve been comparing the properties in the rasterlayer (rasterized vector-layer) with another well working rasterlayer from the North Carolina-datapackage. Question is if I´m converting (rasterizing) the vectorlayer in the wrong way? Here I list the examples of both layers,- 1 from the NC-datapackage (working) and 2, my own file (not working). I see some differences there, but what am i doing wrong?

1, working:

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
C:\Users\Morck\Dropbox\Arbejde\GIS-tutorials\GrassDATABASE\Diverse\ElevationKopi
AREA_OR_POINT=Area
Band 1
COLOR_TABLE_RULE_RGB_0=5.557879e+001 7.572901e+001 0 191 191 0 255 0
COLOR_TABLE_RULE_RGB_1=7.572901e+001 9.587922e+001 0 255 0 255 255 0
COLOR_TABLE_RULE_RGB_2=9.587922e+001 1.160294e+002 255 255 0 255 127 0
COLOR_TABLE_RULE_RGB_3=1.160294e+002 1.361797e+002 255 127 0 191 127 63
COLOR_TABLE_RULE_RGB_4=1.361797e+002 1.563299e+002 191 127 63 20 20 20
COLOR_TABLE_RULES_COUNT=5
STATISTICS_MAXIMUM=128.96199035645
STATISTICS_MEAN=127.18390324034
STATISTICS_MINIMUM=124.20372009277
STATISTICS_STDDEV=1.2340091242347
Dimensions
X: 700 Y: 750 Bands: 1
Origin
638300,220750
Pixel Size
1,-1
No Data Value
nan
Data Type
Float32 - Thirty two bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs
Layer Extent (layer original source projection)
638300.0000000000000000,220000.0000000000000000 : 639000.0000000000000000,220750.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

2, not working

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
C:\Users\Morck\Dropbox\Arbejde\GIS-tutorials\GrassDATABASE\Diverse\VandløbstraceRasterized.tif
AREA_OR_POINT=Area
Band 1
STATISTICS_MAXIMUM=1
STATISTICS_MEAN=0.45934734513274
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=0.49834462137385
Dimensions
X: 226 Y: 256 Bands: 1
Origin
695219,6.18746e+06
Pixel Size
1,-1
No Data Value
*NoDataValue not set*
Data Type
Float64 - Sixty four bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Layer Extent (layer original source projection)
695219.1304354381281883,6187204.6521555976942182 : 695445.1304354381281883,6187460.6521555976942182
Band
Band 1
Band No
1
No Stats
No stats collected yet


texts: Verify if a WMS layer exists or not in GeoServer (OpenLayers3)

texts:
WMSsource = new ol.layer.Image({
    source: new ol.source.ImageWMS({
        url: 'http://localhost:8080/geoserver/alger_data/wms',
        params: {
            'LAYERS': 'alger_data:hotel'
        },
        serverType: 'geoserver',
        crossOrigin: "Anonymous",
        outputFormat: "application/json",
    })
});

map.addLayer(WMSsource);


texts: How to view a single Mapbox Vector Tile?

texts:
I know of setups that display sets of Mapbox Vector Tiles packaged inside a .mbtiles file in a webmap (e.g. http://osm2vectortiles.org/docs/getting-started/ ).

But as part of some experiments I currently have a single Mapbox Vector Tile and need to preview it (as an image) to check its contents.

Does anyone know of a way to do this?


texts: Why ST_Distance give something different to 0?


texts: SDE FeatureClass Granting and Revoking Privilege with ArcObjects

texts:
How can I with ArcObjects GRANT or REVOKE Privileges like SELECT, UPDATE on SDE FeatureClasses to/from Roles?

The equivalent in normal Oracle would be something like this:

GRANT SELECT ON SCHEMA.TABLE TO SOME_ROLE;
REVOKE SELECT ON SCHEMA.TABLE FROM SOME_ROLE;

I can't seem to find the interface to perform this task, and it may not exists?

My dream would be to have an Interface like this:

public interface IPrivilegeManager
{
    void Grant(string privilegeName, string schemaOwnerName, string featureClassName, string roleName);
    void Revoke(string privilegeName, string schemaOwnerName, string featureClassName, string roleName);
}

Related Q: For a SDE FeatureClasses list all Roles that have been Granted any Privileges on it, and which Privileges each Role have on it


texts: Add points to a polygon based on polygon data

texts:
Im working with QGIS 3.x, and i have a polygon layer, and each polygon has a value.
i want create a point layer (like random points) but based on each polygon value.
For example: polygon 1 value = 3 then add 3 random points, polygon 2 value = 12, add 12 random points and so on.
Any lead?


texts: Calculate the values of range, sill, nugget in ordinary kriging?

texts:
I am working on temperature data, and I want to interpolate the data to any location that I specify.

My programming language is asp.net4 C#. I don't know how can I determine the sill, range, nugget for experimental variogram or theoretical variogram.

I make the following code to calculate these variables in experimental variogram:

public void Calculate_EmpericalVari_nugget_sill_range()
{
double H=0,VaiogramValue;
int N = KnownPoints.Count;
for (int j = 0; j < N; j++)//calculate the emperical variogram
{
for (int i = 0; i < N; i++)
{
H = AucledianDistance(KnownPoints[i].X, KnownPoints[i].Y,
KnownPoints[j].X,     KnownPoints[j].Y);
VaiogramValue = 0.5 * Math.Pow((KnownPoints[i].Z - KnownPoints[j].Z), 2);
Emper_Variogram[i][j] = VaiogramValue;
if (H==0) nugget = VaiogramValue;
if (sill < VaiogramValue)
{
sill = VaiogramValue;
range = H;
}
}
}
}

I work according to the following that sill is the maximum value and range is the distance when I get the sill and nugget is the value at distance equal zero?

Is this right? And do they remain constants when I calculate the theoretical semivariogram?


texts: How to use Expressions in QML file in QGIS 3.x

texts:
I have a plugin that uses optional QML files for layer labels.  If I use a standard field name for the value, it works fine, but if I use an expression, it will not display the label UNLESS I open and close the Properties dialog for the layer.
This works:
<property key="labeling/fieldName" value="id"/>
This does not:
<property key="labeling/fieldName" value="concat(id, ' - ', acquisition)"/>
The code is fairly simple:
        style_path = os.path.normpath(QgsApplication.qgisSettingsDirPath() +
                                      f'/python/plugins/S2SRoadNet/{fnm}')
        self.dbshp.segments_layer.loadNamedStyle(style_path)
        self.iface.layerTreeView().refreshLayerSymbology(self.dbshp.segments_layer.id())

What is the Properties dialog doing that refreshLayerSymbology is not and how can I run whatever that process is in code to get the expression to work in the QML file?


texts: How to download whole ECW or JPG tiles from web viewer?

texts:
I've found an online server which has some aerial tiles I'd like to save at maximal resolution for offline use, I'd even be happy with just the whole JPG picture to be later georeferenced again by myself.

I did try to test the source URL and found it's probably using an ERDAS ecw server setup and JPG tiles (also probably found the remote ECW directory but seemingly gives an access error message), however being no more experienced than reading some HTML, I don't know how I'd download the whole image at maximum zoom.

Upon looking potentially similar discussions I came across a software called taho, however I think it may not apply to my case.

Could you help me?

Here is the webpage viewer URL: //fotogrammi.regione.abruzzo.it/sampleiws/fotogrammi.asp?VOLO=VIGM54&&LOTTO=1&&STR_=48&&FOTO=1686

and here a sample tile URL //fotogrammi.regione.abruzzo.it/ecwp/ImageX.dll?image?cache=true&transparent=true&type=jpg&l=0&tx=9&ty=9&ts=256&fill=FFFFFF&quality=60&layers=/VIGM54/1/48/1686.ecw&srs=EPSG:RAW

which, also considering the peculiar "error" page (shows a local "C:\Program Files\ERDAS\Image Web Server" remote directory line) that loads up, led me to think that the ecw is stored at .../ecwp/VIGM54/1/48/1686.ecw


texts: PostGIS layer displays correctly in QGIS, generates error in ArcMap?

texts:
I have a road centerline layer in PostGIS that displays correctly in QGIS 2.18.10, but errors in ArcMap 10.5.1 with

  The specified coordinate exceeds the valid coordinate range

Other PostGIS layers show correctly in ArcMap.

What is wrong with this road centerline layer that ArcMap doesn't like but QGIS doesn't mind?


texts: Deleting values from raster file using ArcGIS Raster Calculator?

texts:
I have a raster that denotes certain landcover types.

How do I delete all of the pixels except for the ones I want?

I tried going through the symbology tab and while that deletes it for display it doesn't change the raster file itself, so if I were to delete and reload the file it would revert back.


texts:
I want to make interactive maps online/on a server but not sure whether leaflet or MapGuide program is better for me to delve into.

Which software will allow me to render resource intensive layers quickly?
Which software is most compatible with QGIS?
Which software would allow users to switch between up to 30 mutually exclusive layers?
Which software is easier to update a couple of layers without having to re-do to much?


texts: GeoTools Oracle spatial plugin samples / documentation

texts: I plan to use GeoTools and I saw that there is an Oracle plugin available, but the documentation is not very detailed.  Are there any examples or documentation available to access the functionality of the plugin?


texts:
I've started getting errors from geoserver when trying to request WMS layers.
I'm running it on windows server (the OpenGeo suite version)

The errors are :

WARN [mortbay.log] - handle failed
java.io.IOException: FULL head
    at org.mortbay.jetty.HttpParser.parseNext(HttpParser.java:276)
    at org.mortbay.jetty.HttpParser.parseAvailable(HttpParser.java:205)
    at org.mortbay.jetty.HttpConnection.handle(HttpConnection.java:380)
    at org.mortbay.io.nio.SelectChannelEndPoint.run(SelectChannelEndPoint.java:395)
    at org.mortbay.thread.BoundedThreadPool$PoolThread.run(BoundedThreadPool.java:450)

And I'm getting pink tiles on FF or empty tiles on Chrome/IE.

The weird thing is that on geoserver itlsef I can preview the layer.

any ideas?


texts: Align existing points along a line with even spacing

texts:
I have a number of existing points which are snapped to a line, however they are unevenly spaced along the line, like this:

How do I align them so they are evenly spaced out along the line like this:

I know how to create new, evenly spaced out points along a line, however these are existing points which have attribute info that I want to keep.
I am using QGIS 3.14.16.


texts: How can I split a polygon into its sub polygons?

texts:
I am trying to define user generated polygons with Esri's ArcGIS Engine 10.1.
Here's how it works at the moment:

User clicks on the map and creates points
User connects the points

I have a collection of line segments (IPointCollection) in memory.  The map looks something like this.

I want to input all my line segments into an object.  Then I need the object to tell me that I have 2 polygons, and one of the polygons has 5 sub polygons.

I created a ISegmentCollection with all the segments in both polygons, then I cast the  ExteriorRingBag into an IGeometryCollection, and I have my 2 main polygons.

The bounty goes to whoever can tell me how to take the IPolygon on the left, and get it to return it's 5 sub polygons.


texts: Translating "Choose Multiple" Checkbox Forms from Survey123 to ArcGIS Online Maps?

texts:
I have a request for a few ArcGIS Online forms that require possible multiple answers in certain situations. I did notice that this is can be done using Survey123 forms and checkboxes.

But when the same feature is being edited in an ArcGIS Online Map it does not show the checkboxes, it only shows a large textbox. I was curious if there is a workaround that allows you to be able to choose mutiple "answers" in arcGIS Online similar to how it looks in Survey123. There are some people that like using Survey123, but a majority of people prefer using Collector or the ArcGIS Online webmaps to enter the data. Here are a couple images to show what I am referring to.

Form in Survey123

Form when Seen in ArcGIS Online


texts: Retrieve Google Maps API friendly latitudes and longitudes from shapefile polygon

texts:
I have a shapefile with a polygon in it, and I am trying to the latitudes and longitudes associated with it so I can use them with the Google Maps API.

So far I have been working with QGIS Lyon and have managed to get some sort of coordinates, but they are not compatible with Google Maps. This is what the coordinates look like in the CSV:

"shapeid","x","y"
"0.0","543905.4","183199.1"
"0.0","543905.0","183198.9"
...

I have used a plugin called MMQGIS to get the values in the CSV file.

I have also tried to use Vector->OpenStreetMap->Download Data, and this does give coordinates in the right format, but the output file is very large, 12MB, and it has around 92k nodes, which is definitely way too many.

I've come to a point where I can't think of anything else and Google isn't being very helpful either (or maybe I'm not doing the right search).

So, it would be very helpful if someone could point me to the right direction.


texts: How to estimate flow paths between contaminant sources and catchment offtake (exposure point) at catchment?

texts:
I am trying to estimate flow paths between contaminant sources and the catchment offtake (exposure point) at the catchment.

Can anybody give me some help in how to proceed?


texts: Where does the Style Manager store the files in and Why does the Style Manager does not show some edit color ramps?

texts:
Where are Color Ramp files saved? Every time I make edits, it does not show up on the Style Manager. Does it save it to the hard drive or in the Style Manager?
I need more information.

Thanks


texts:
I have loaded data from a MapInfo TAB file into PostGIS using EasyLoader. I can view the data in MapInfo now, except that the feature-based styles are being ignored.

When I look at the map catalogue, I can see that it specifies the MI_STYLE column  as the column to use for styling, but this seems to be ignored by MapInfo.

I'm using MapInfo v15, Postgres 9.5 and PostGIS 2.2.


texts:
Is it possible to get search results from ArcGIS Online and ArcGIS OpenData as JSON?

When I search for, say 'population', the URL generated is either http://www.arcgis.com/home/search.html?q=population&t=content or http://opendata.arcgis.com/datasets?q=population

Can I append something to the URL to get a JSON response rather than an HTML page?


texts: VB scripting in ArcGIS Field Calculator?

texts:
In my attribute tables field I am trying to get rid of the "&" in between names i.e. "John Adams & wife" and want it to say "John Adams and Wife" in all the attributes with "&".

What is the VB for that?

I can't seem to find it in the help page of ArcGIS.


texts: How to only select lines (geometry features) going through a polygon?

texts:
I'm currently trying to select bridges (line features), from a shapefile containing roads, which go across water (a polygon). The problem is when i try to select features crossing or intersecting the water polygon (both in arcgis and qgis), it includes roads on piers as well. It is not possible to manually select the features, as the dataset consists of approximately 1.5 million features.

What i have noticed is the lines on piers simply intersect/crosses the water polygons, while bridges i have found goes "through" the polygon (one part of a line feature enters through one polygon boundary, but the line feature also exits either through the same boundary or another boundary).

Is there any spatial operator/plugin that can select features going through the polygon instead of intersecting/crossing the polygon? basically a double-crossing spatial operation.

edit: the road dataset does not contain any attributes to distinguish between different road "types"
edit: As richard pointed out, it probably needs some graphic examples. Below are some simple visualizations of the red lines only crossing the green water polygon once (left) and red lines crossing the water polygon twice. I'm trying to get the double-crossing features.


texts:
I am trying to create a polygon showing crown spread of a tree. Using QGIS.

I am aware you can achieve this with geometry generators, but I need to display symbols in legend and looking for a polygon creator so spatial analysis can be carried out in the future as well.

Does anyone know how to create a crown spread polygon or convert the geometry generator to a polygon?


texts: KML created by Google Earth Engine can't be read on Google MyMaps

texts:
I have a FeatureCollection in Google Earth Engine, and I'm exporting it to a kml (and storing on Google Drive) using the code below:

// Export the FeatureCollection to a KML file.
Export.table.toDrive({
  collection: fcPR,
  description:'citiesOfParanaToDrive',
  fileFormat: 'KML'
});

The permissions are set to Public (anyone can read and access the data).

And when I try to import it in Google MyMaps, I get the message: "The file contains invalid or not supported data, or is bigger than supported". Opening in Google Earth (installed application on computer) works.

Reading the documentation, the max number of lines is 2000, but my file contains only 400.

Is there some step that I am missing?


texts: How to render the Sea depth level in different color?

texts:
I would Like to adapt an old code in c++ using the api qgis 1.6 :

const QgsRenderer* lrenderer = pLayer->renderer();
if( lrenderer == NULL )
    lrenderer = new QgsSingleSymbolRenderer( pLayer->geometryType() );
QList<QgsSymbol*> lSymbols      = lrenderer->symbols();
QgsSymbol* lsy = lSymbols.at(0);
lrenderer = new QgsUniqueValueRenderer( QGis::Polygon );
QStringList lListValue; lListValue  << "0" << "3" << "5"<< "10";
for(int i = 0 ; i < lListValue.size() ; i++) {
    QgsSymbol* newsymbol = new QgsSymbol( *lsy);
    newsymbol->setLowerValue( QString::number(i) );
    newsymbol->setFillColor(QColor(115, 182, 239));
    newsymbol->setColor( QColor(125, 137, 140) );
    newsymbol->setLineWidth(0.3);
    ((QgsUniqueValueRenderer*)lrenderer)->insertValue( lListValue.at(i), newsymbol ); //A
}
((QgsUniqueValueRenderer*)lrenderer)->updateSymbolAttributes(); //B
((QgsUniqueValueRenderer*)lrenderer)->setClassificationField(lMapField.value("DRVAL2").toInt() ); //C

This is a qgis api 1.6 c++ code
Now I want to do the same thing using the qgis 2.8.2 api c++
I have started to code it :

lrenderer = QgsFeatureRendererV2::defaultRenderer( QGis::Polygon );
QList<QgsSymbolV2*>         lSymbols    = lrenderer->symbols();
QList<QgsSymbolV2*>         lSymbolsToRender    = lrenderer->symbols();
QStringList lListValue;
lListValue  << "0" << "3" << "5"<< "10";

for(int i = 0 ; i < lListValue.size() ; i++) {
    QgsSymbolV2* newsymbol = lSymbols.at(0);
    QgsSymbolLayerV2List  lSymbolLayerList = newsymbol->symbolLayers();
    QgsLineSymbolV2 *MyLine = new QgsLineSymbolV2(lSymbolLayerList);
    MyLine->setWidth(0.3);
    QgsSymbolLayerV2 * lSymbolLayer;
    lSymbolLayer = lSymbolLayerList.at(0);
    lSymbolLayer->setFillColor(QColor(115, 182, 239));
    lSymbolLayer->setColor( QColor(125, 137, 140) );
    newsymbol->insertSymbolLayer(lListValue.at(i).toInt(), lSymbolLayer);
    lSymbolsToRender.append(newsymbol);
}

I don't know how i can code the A, B and C comment of the old code ?
Does anyone have sample on how to renderer the deep of the sea with different color, ie level 0 is light blue 10 feet deep is blue and more is dark blue


texts: QGIS: Using OSgeo4W shell to login to and create a new PostGIS database

texts:
I'm trying to create a new database through the OSgeo4W shell, as opposed to loading QGIS and hitting the 'Manage DB' menu to make one an pipe in SQL scripts. I haven't seen any documentation on this specific issue.
I've tried piping in psql commands with the -p , -U and createdb commands to login and attempt to make a new database. I haven't had any success so far. I've mainly gotten connection refused messages, even though I type in the login information.
I'm not sure what the issue is here, or what I may be forgetting to do.
Shown are a couple attempts:


texts: Query to Calculating Rail Length Using Overpass API / OSMnx Python

texts:
I'm new to overpass queries. I'm trying to write a python script to pull data on the length of rail in major global cities. I am having a problem in the length of track (using osmnx.basic_stats(G)) seems abnormally large.

I have a two questions:

For a line of track (running in two directions side by side), it looks like overpass returns two tracks for each direction (eg, four lines of track in total). See this screenshot below in London.

However, for a subway, overpass returns a single line of track in each direction. See screenshot below.

To calculate the total length of track then, is the best way to do this simply to collect these differently and account for these differences by dividing the length of rail track by 2?

My second question is, what would be an appropriate query to remove disused tracks and industrial tracks. You can see the OSM Railway tag guide here.

FYI, my code to collect the data is as follows:

import osmnx as ox

# Pull multigraph using Londons lat/long
G = ox.graph_from_point((51.5073509,-0.1277583),
                        distance = 20000,
                        distance_type = 'bbox',
                        infrastructure = 'way["railway"~"rail"]',
                        network_type = 'none',
                        )

# Print length of track in KM
print(ox.basic_stats(G)['edge_length_total'] / 1000)

# Visualise the data that has been pulled
ox.plot_graph_folium(G)


texts: OpenLayers 5 and GeoJSON with two geometry attributes per feature

texts:
I have a GeoJSON file in which each feature has two geometry attributes, one is the default geometry property defined by the GeoJSON standar, and the other one is a property under the standarized properties entry. My intention is, render each geometry with a different style.

Following you can find an example of a GeoJSON feature, and the attempt to define the OpenLayers style.

{
"type": "FeatureCollection",
"crs": {
    "type": "name",
    "properties": {
        "name": "CRS:84"
    }
},
"features": [
    {
        "type": "Feature",
        "geometry": {
            "type": "Polygon",
            "coordinates": [
                [
                    [
                        -13.7331,
                        43.0833
                    ],
                    [
                        -12.6167,
                        41.85
                    ],
                    [
                        -19.7833,
                        37.8833
                    ],
                    [
                        -20.9333,
                        39.0333
                    ],
                    [
                        -13.7331,
                        43.0833
                    ]
                ]
            ]
        },
        "properties": {
            "type": "StaticAre",
            "name": "nameValue",
            "test-attr1": "test-attr1",
            "geom_alt": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [
                            -13.7333,
                            43.0833
                        ],
                        [
                            -12.6167,
                            41.85
                        ],
                        [
                            -19.7833,
                            37.8833
                        ],
                        [
                            -20.9333,
                            39.0333
                        ],
                        [
                            -13.7333,
                            43.0833
                        ]
                    ]
                ]
            }
        },
        "id": "6764"
    }
]

}

var style = [
new ol.style.Style({
    stroke: new ol.style.Stroke({
        color: 'red',
        width: 2
    })
}),
new ol.style.Style({
    geometry: 'geom_alt',
    stroke: new ol.style.Stroke({
        color: 'green',
        width: 2
    })
})
];

I think that the problem is related with that openlayers does not identified the 'geom_alt' property as a geometry.

Do you know if there is a way of doing this?


texts: When PyQGIS script is run through console layer does not load successfully

texts:
I made a small script to save a shapefile to disk. I want to run the script through console. When I run the script through QGIS IDE, it runs successfully. But when I run it through the QGIS Python command prompt it says layer failed to load successfully.

What am I doing wrong?

from qgis.core import *
import qgis.utils

_layer = QgsVectorLayer("C:/Arundhati/RadioSignal/MasterLine.shp", "MASTERLINE", "ogr")

if not _layer.isValid():
  print("Layer failed to load!")
else:
    print("Layer was loaded successfully!")

_writer = QgsVectorFileWriter.writeAsVectorFormat(_layer,
                               "C:/Arundhati/RadioSignal/MasterLine12.shp",
                               'utf-8',
                               QgsCoordinateReferenceSystem(),
                               'ESRI Shapefile')

I run the above file through QGIS Python console like the following:

Additionally Python also crashes:

What Am I doing wrong?

How to run the same script from outside the QGIS environment.


texts: Capture click while drawing a LineString in openlayers and query GeoServer WMS

texts:
While drawing a LineString in OpenLayers version 4.6.5, I would like to hook into the clicking of each point, and make an ajax call to a backend WMS service (geoserver).

If a Feature on the Tile layer is not clicked, reject the plotting of the point - like connect the dots of the features in the TILE layer.

I already have an asynchronous method that will retrieve the feature based on coordinates from the WMS service, returning a Promise of the Feature.

Here's some psuedo-code to explain what I want to do.  I think it will be more difficult because it is asynchronous and requires a callback.

function onPointCreate(obj) {
  getFeatureFromWMS(obj.coordinate).then(
    function(feature) {
      if(feature) {
          // allow plotting of the point in the polyline
      } else {
          alert("Sorry, you must click on a feature");
          // don't allow plotting
      }
    });
}

I just can't find where to hook something like this up. I've tried geometryFunction and the drawstart event without success.


texts:
I'd like to add some Strava WMTS to ArcGIS.

The QuickMapServices QGIS plugin can achieve this.

The layer details of this plugin show this information with the following:

URL: http://globalheat.strava.com/tiles/cycling/color1/{z}/{x}/{y}.png
Z min: None
Z max: None
Y origin top: True

Is there any way to use the information given to add this service in ArcMap?


texts: How to keep the cummulative polygons in time manager for QGIS?

texts: I use time manager for a polygon layers with dates and it works. I would like to see the cummulative addition of my polygons. So I have 4 periods I want the poylgons from 1900-1950 to be added to the new ones in 1950-2000 etc in the next time step. How can I perform that or does another open source tool exist for that?


texts: Calculating distance between points from CSV-file in QGIS

texts:
I have a CSV-file with latitudes and longitudes and I should calculate the distance between the points in ascending order (distance from 1 to 2, from 2 to 3, etc.)
I tried with the distance matrix but it calculates each point the distance with the others.
In addition, I am working in EPSG:23033 - ED50 / UTM zone 33N but when I made the distance matrix no measurement is returned. The ruler tool tells me that EPSG:7022 works automatically but I do not want to work with that CRS.
My initial CSV-file with 241rows looks like this:

I would like a new column.


texts: Order of polygon vertices in general GIS: clockwise or counterclockwise

texts:
Two days ago I asked a question about the internal storing order for the vertices of a polygon in ESRI shapefiles. That question was answered (Are polygons stored clockwise or counterclockwise in a shapefile?) and it was also answered in an old post (Polygon creation (Clockwise rotation or not))
But now my question is more general, and I don't know if it has an unique answer.
Is the clockwise order only for ESRI shapefiles or for general GIS formats?
And what about the internal representation for a GIS software?
For example, if I use QGIS and I read a *.shp containing polygons I suppose the internal representation of the outer bound is clockwise as in the original shapefile, but what about for all the file formats supported by QGIS?
And for ArcGIS?
And in the case there exists a file format with polygons stored in counterclockwise, if these files are loaded in QGIS, ArcGIS, etc., is the orientation changed internally, so if I read the data using PyQGIS, for example, the polygons are clockwise ordered?
My purpose is to write a plugin for QGIS, but the source of data can be ESRI shapefiles or other formats. As I need to check the angles between consecutive sides of polygons using their azimuths I need to known if the order is clockwise. One solution is computing the area of each polygon and, if I remember correctly, if it is positive the order is clockwise and if negative the order is counterclockwise.
Area computation is not an intensive task, so it will not slow my plugin so much.
But in the special case of QGIS, anybody knows if it stores the polygons clockwise or counterclockwise, regardless the order in the original source?
By now I'm working with ESRI shapefiles and the coordinates in layer.getFeatures().geometry().asPolygon() are stored in clockwise for the outer border and counterclockwise for the holes, i. e. as in the original *.shp.


texts: Database View in Oracle 9i with SDE 9.3.1

texts:
If I have a sql statement creating a view against ArcSDE 9.3.1 on Oracle 9i, is the best practice to create the view in ArcCatalog just creating the name (kind of shell of the view) then going in to PLSQL to modify the view?  What are the best practices on this?

Reading this.. http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=Using_database_views ..I am still left without great advice.  Thanks.

Here is the view

CREATE OR REPLACE VIEW SERIES_DETAIL AS
SELECT LineLoop.LineName, Cont_StationSeries.EventID AS RouteID, Eng_StationSeries.EventID AS SeriesID, Eng_StationSeries.SeriesOrder, Eng_StationSeries.BeginStation, Eng_StationSeries.EndStation, Begin_Cont_ControlPoint.StationValue AS BeginMeasure, End_Cont_ControlPoint.StationValue AS EndMeasure

FROM (((((LINELOOP LineLoop INNER JOIN STATIONSERIES Cont_StationSeries ON LineLoop.EventID=Cont_StationSeries.LineLoopEventID) INNER JOIN STATIONSERIES Eng_StationSeries ON LineLoop.EventID=Eng_StationSeries.LineLoopEventID) INNER JOIN CONTROLPOINT Begin_Eng_ControlPoint ON (Eng_StationSeries.BeginStation=Begin_Eng_ControlPoint.StationValue) AND (Eng_StationSeries.EventID=Begin_Eng_ControlPoint.StationSeriesEventID)) INNER JOIN CONTROLPOINT End_Eng_ControlPoint ON (Eng_StationSeries.EndStation=End_Eng_ControlPoint.StationValue) AND (Eng_StationSeries.EventID=End_Eng_ControlPoint.StationSeriesEventID)) INNER JOIN CONTROLPOINT End_Cont_ControlPoint ON End_Eng_ControlPoint.GroupEventID=End_Cont_ControlPoint.GroupEventID) INNER JOIN CONTROLPOINT Begin_Cont_ControlPoint ON Begin_Eng_ControlPoint.GroupEventID=Begin_Cont_ControlPoint.GroupEventID

WHERE (((Cont_StationSeries.SubTypeCD)=1) AND ((Eng_StationSeries.SubTypeCD)=2) AND ((Begin_Eng_ControlPoint.SubTypeCD)=2) AND ((End_Eng_ControlPoint.SubTypeCD)=2) AND ((Begin_Cont_ControlPoint.SubTypeCD)=1) AND ((End_Cont_ControlPoint.SubTypeCD)=1))

ORDER BY LineLoop.LineName, Cont_StationSeries.EventID, Eng_StationSeries.SeriesOrder


texts: I'm a beginner QGIS user and am struggling with creating a grid. I've used the MMQGIS geometry import function to import data from the CSV file shown below but they pop up as points on QGIS. I want to make it into a grid so that I can pretty much have each of the points be in a box format. I've also tried to use the join function but am not sure what variable both the grid and the points have so that I can join them.


texts: Fastest way to union a set of polygons in R

texts:
I want to find the fastest way to union a set of polygons into one large polygon.

Lets first get some data:

# Load libraries
library('raster')
library('geosphere')
library('mapview')
library(maptools)
library(rgeos)
library(sf)

# Get SpatialPolygonsDataFrame object example
pols<- getData('GADM', country = 'DK', level = 2)

#Project to suitable projection (to be able to calculate area, see later
utm32 = "+proj=utm +zone=32 +ellps=WGS84 +units=m +no_defs"
pols<- spTransform(pols, CRS(utm32))

mapview(pols)

# 1st approach: maptools::unionSpatialPolygons
system.time(pol1 <- unionSpatialPolygons(pols,rep(1, length(pols))))
#  bruger   system forløbet
#    3.67     0.03     3.72

# 2nd approach: rgeos::gUnion
system.time(pol2 <- gUnaryUnion(pols, id = pols@data$NAME_0))
#  bruger   system forløbet
#    3.69     0.00     3.74

#3rd appraoch: sf:st_union
pols_sf <- st_as_sf(pols)
system.time(pol3 <- st_union(pols_sf))
#  bruger   system forløbet
#    3.67     0.02     3.68

# 4th approach: rgeos::gBuffer
system.time(pol4 <- gBuffer(pols, byid=F, width=0))
#  bruger   system forløbet
#    1.13     0.00     1.16

Of the four approaches, the three first is very similar, whereas #4 is significantly faster. My problem is that the polygons are not identical:

identical(pol1, pol4)
[1] FALSE

And the areas are slightly different:

paste(area(pol1))
[1] "43122105144.9307"

paste(area(pol2))
[1] "43122105144.9307"

pol3 <- as(pol3, "Spatial")
paste(area(pol3))
[1] "43122105144.9724"

paste(area(pol4))
[1] "43122105144.9062"

Why is this, and is there a reason for using one approach over the other (apart from processing time)?
Also, do you know of any approaches that are faster?

EDIT:

I did some more testing with more polygons, and it seems as method 1-3 only gets slightly slower with larger dataset, whereas method 4 gets very slow.


texts: Attribute field names truncation when creating shapefiles

texts: I am wondering what is the expression that truncate a given name to the 10 characters name used in ESRI shapefiles when an object is saved as shapefile with the function writeOGR in R. For example, the attribute field from a spatial dataframe "Ablepharus kitaibelii" is converted to "Ablphrk" when saved as shapefile, how the truncation is done?


texts: ArcGIS Desktop can't return concurrent license

texts:
We have one concurrent use ArcView license, managed by a license server box in the office. The server reports that one of our desktops has the license 'borrowed', but when looking at the Desktop Admin tool on that box, it doesn't show anything about having a license checked out. It seems the license is somewhere lost in the ether.

Is there a way, from the desktop side, to 'force return' or destroy the local license file? There appears to be no way for me to get the license returned from the desktop side. I've filed a report with Esri, but not holding out hope that my issue will be fixed ASAP.

Anyone here have experience with this problem? Would a reinstall of the desktop software help?


texts: QGIS 'Export layer as' retaining associated joins and widgets

texts:
I have a layer with multiple joins and widgets (value map). This is a test layer and I would like to export it so that it can be renamed multiple times and used in different projects however when I export it it loses all the joins and widgets. Is there a solution for this?

Thanks!


texts: Pure PHP Library for Geometry Operations 7

texts: I am working on a php project on a shared webhost. so I am limited into mysql and php. I want to do some geometry operations in php side. In fact it is a routing algorithm which works in php. I need a library to union and intersect geometries using php. I know about GeoPHP but it needs GEOS and it is not possible to install it on shared hosts. Is there any geospatial php libraries?


texts:
I need to remove duplicates from a field named "Intersecti"

Here is my code.  I used list..set. It's not removing any duplicates.

import arcpy
duplicates = "G:\\xStreetNew\\Duplicates.shp"
duplicates_List = []
with arcpy.da.SearchCursor(duplicates, ['Intersecti']) as Cur:
 for row in Cur:
    duplicates_List.insert(0,row[0])
    duplicates_List = list(set(duplicates_List))
 print duplicates_List


texts: What's the name of the area between two geographic coordinates?

texts: Is there a name for the shape of the area defined by two sets of geographic coordinates, for example by a northwestern and southeastern corner? In other words, the shape bounded by two degrees of longitudes and two degrees of latitude? Since it's on a sphere/ellipsoid, it's not a rectangle or a trapezoid. Is there another good name for it then? A quadrangle? I've looked on wikipedia, google and stackexchange, but couldn't find anything.


texts: MapServer Runtime substitution for geotransform

texts:
I am trying for a way to get a variable geo transformation of a few line strings from a shapefile:

VALIDATION
  'TIME' '^[0-9\-]*T[0-9:]*\.0Z$'
  'SMOOTH' '^[0-9]*$'
  'SMOOTHITER' '^[0-9]*$'
END
GEOMTRANSFORM (smoothsia([shape], %SMOOTH%, %SMOOTHITER%))
LABELITEM "value"
CLASS
  EXPRESSION ([value] <= 100)
  STYLE
    WIDTH 2
    COLOR 203 35 5
  END
  LABEL
    PRIORITY 1
    PARTIALS FALSE
    ...

As I finally found in the documentation this does not seem possible as I want. I would be fine with defining the layers several times but then I need to copy the whole color table and all the info over and over and over in the mapfile. So it would be nice to get this done with a parameter a user can give to the WMS.

Is there an easier way then copy paste of the full layer just changing the two parameters? Also this would make a hard coded selection which I would like to avoid.


texts: Hansen Dataset - Extracting Yearly Forest Cover from Hansen

texts: Could someone give an suggestion about, how to extract yearly forest cover from Hansen Dataset? E.g forest cover 2000, forest cover 2001, and so on.


texts: Identifying pixels with sun glint in a satellite image (for atmospheric correction)

texts: Wondering if anyone knows a way (perhaps objective) to identify pixels that are infected with sun glint in an image. I'm using Hedley's 2005 approach in a python to deglint my image but this approach requires selecting pixels over deep waters with glint/specular reflection. I'm having trouble identifying this visually in ENVI - even using different contrast options.


texts: Calculating polygon area within polygon using ArcGIS Pro

texts:
I have a soil layer overlaying a watershed layer in ArcGIS Pro.
How do I calculate the different types of soil area within the watershed?
I have tried summarized within, but it does not work.


texts:
I use 9 base map layers from .SHP, .SLD, .DBF and etc. files. Their common HDD-size is 134 MB. And I faced with a perfomance issue on my PC. It takes aboout 20-25 seconds to zoom or to move the map, App holds for this period and it's not good at all. Also I want to note that the larger the scale, the more the App slows down. I think it's because it has to calculate more different elements and features and when we zoom their number reduces. At the city-building zoom level it works pretty fast.

So is there some tricks or approaches for performance improvments?

Now for temporal solution I hide\show 8 of 9 layers and it reacts quicker but all the same it takes the same time to redraw the map. Of course I use the parallel thread for my draw method to avoid the other App's GUI and scenes holds. For now I didn't try to use the tiles funcs yet.

May be tiles will be some kind of solution or maybe there are other
aprroaches?
Can you suggest articles or code examples for this topic?

Except these ones:
- doc1
- doc2

Here is my code

This method draws the map:

public void drawMap(GraphicsContext gc) {
    if (repaint) {
        repaint = false;
        StreamingRenderer draw = new StreamingRenderer();
        draw.setMapContent(map);
        FXGraphics2D graphics = new FXGraphics2D(gc);
        graphics.setBackground(colorBackground);
        graphics.clearRect(0, 0, (int) canvas.getWidth(), (int) canvas.getHeight());
        Rectangle rectangle = new Rectangle((int) canvas.getWidth(), (int) canvas.getHeight());

        draw.paint(graphics, rectangle, map.getViewport().getBounds());
    }
    else {return;}
}

This method sets new display area:

public void doSetDisplayArea(ReferencedEnvelope envelope) {
    map.getViewport().setBounds(envelope);
    repaint = true;
}

Here is my hide method's code:

public void hideLayers(int indexFrom) {
    //Hide layers for performance achievments
    java.util.List<Layer> layers;
    layers = map.layers();
    for (indexFrom = indexFrom; indexFrom < layers.size(); indexFrom++){
        layers.get(indexFrom).setVisible(false);
        Layer l = layers.get(indexFrom);
        System.out.println("SYSTEM: layer " + l.toString() + " was hidden. ");
    }
}

and here is the show method's code:

public void showLayers(int indexFrom, ReferencedEnvelope envelope) {
    //Show layers for performance achievments
    java.util.List<Layer> layers;
    layers = map.layers();
    for (indexFrom = indexFrom; indexFrom < layers.size(); indexFrom++){
        layers.get(indexFrom).setVisible(true);
        doSetDisplayArea(envelope);
        Layer l = layers.get(indexFrom);
        System.out.println("SYSTEM: layer " + l.toString() + " was shown. ");
    }
}


texts: How to persist a polygon to SQL Server database and vice versa?

texts:
How do I extract all the points from a Polygon drawn in ArcMAP and then persist that data into a SQLServer database table containing a geometry column and then when I've got the data persisted read from the database and draw the polygon again.

I'm working with ArcGIS 10 sp1 using ArcObjects for .NET (C#) sp1. I've made an add-in tool to draw the polygon and want that add-in to persist it too. I know how to do the last steps e.g putting the geometry data into the database and drawing an IPolygon instance in ArcMAP, what I don't know how to do is to somehow "convert" the Polygon into something that will "fit" into the geometry column of the database and vice versa.

I'm completely new to ArcObjects programming and GIS systems in general, started two weeks ago, so bare with me...


texts: Error when creating enterprise geodatabase

texts:
I'm trying to create an enterprise geodatabase.  My problem is when I try to create a corporate geodatabase.  The error is:

  Cannot connect to database because the database client software failed
  to load. Be sure the database client software is installed and
  configured Correctly. Failed to execute (CreateEnterpriseGeodatabase).

Here is what I have tried so far.

I have installed on my server:

Windows Server 2012 R2 64-bit
ArcGIS for Server 10.2.2
ArcGIS Desktop 10.2.2
PostgreSQL 9.2.2

I have also downloaded the following client libraries:

PostgreSQL 9.2.2 64-bit
iconv.dll libeay32.dll libintl-8.dll libpq.dll ssleay32.dll

I have pasted these libraries into the following folders on my server:

C: \ Program Files \ ArcGIS \ Server \ bin
C: \ Program Files (x86) \ ArcGIS \ Desktop10.2 \ bin

I have also pasted the ST_Geometry libraries in the lib directory of PostgreSQL and have edited the text file, pg_hba file.


texts:
I am trying to build the qwc2 demo app from https://github.com/qgis/qwc2-demo-app on a windows machine. I have a working QGIS server and a customized "themesConfig.json".
After "git clone" I do "yarn install" like the documentation suggests, but I get some deprecation warnings:

Do I have to use a specific node or yarn version to compile qwc2?
If I ignore the warnings and continue with "yarn start" the process stops with an error "Missing package 'UUID'":

I tried to add the package "uuid": "8.3.2" to \qwc2\package.json, did "yarn install" and "yarn start" but then it complained about another missing package: "glob".
If I add this package, there are some other missing packages until it complains about missing "chokidar". At this point adding it to packages.json doesnt work and I am stuck.
I am new to npm/yarn.
What am i doing wrong here?
Do i need some special switch to let yarn pull all required dependencies on its own?
Or is there a problem in the qwc2 repository?


texts: Which is the correct path for a hyperlink in an ArcMap Table?

texts:
I am having problems making the hyperlink tool open a hyperlink with a relative path in ArcMap 10.7.
The relative path to the pdf file is in a table joined by a "relate" to a polygon feature class in a FGDB. The Hyperlink setting is set to script with the following code: "PDF" is the field name.
Function OpenLink ( [PDF] )
 Dim path
 Dim objShell
 path = [PDF]
 Set objShell = CreateObject("Shell.Application")
 objShell.ShellExecute path, "", "", "open", 1
End Function

I can successfully open the hyperlink in the info window using the yellow hyperlink tool when the PDF file is saved in the same directory as the MXD but when the PDF file is moved to a directory i.e ../PDF/document.pdf ArcMap does not recognise it as a hyperlink.
I have tried all the following formats with and without quotation marks in the table.
   ../PDF/document.pdf
   /PDF/document.pdf
   PDF/document.pdf

What am I doing wrong?


texts: Performing nighttime light intercalibration of DMSP/OLS in GRASS GIS?

texts:
I have 1 nightlight image of DMSP/OLS. It is uncalibrated. I want to calibrate the image based on Wu's (2013) methodology. Every time I try to run the command "i.nightlights.intercalibration" on GRASS GIS (v. 7.4) I get that error "no data base element found".

Why is that?

I have made every step of Wu's methodology.


texts: MapInfo 'Info' tool - querying by polygon rather than at a point

texts:
MapInfo Professional (12.5) has a tool called 'Info' which interrogates the data for all visible and selectable tables when the user clicks at a specific point. It lists all objects from all tables, and then allows the user to click on each one to view the attributes and values of that object.

Are any MapInfo Pro users aware of a tool which can do this over an area?

For example, showing all objects from all visible tables over a user-defined polygon, or within a user-specified radius from a point?


texts: Getting list of all embassies using Overpass API?

texts:
I'm new to OSM and the overpass API. I've tried using

[out:csv("name:en")];relation["amenity"="embassy"];out;

In order to get the list of all of the embassies mentioned, but the result of the query only produces 38 embassies. That can't be right. Is there something wrong with my query or this is a limitation of the database or something else?


texts: Automating MODIS image download and its processing in ArcGIS Desktop?

texts:
For MODIS satellite images how to automate the download, unzip, import and process to point and then to Excel for each climate satellite MODIS layer.

There are 3650 days of imagery to process with 8 layers per day across modis aqua and terra, giving LST AM and PM and the corresponding timestamps.

So this is 29,200 layers to be run over a set of point files to give the corresponding satellite data at each point for each layer type needed.

To me it seems like a massive undertaking but my knowledge of coding and automation is limited

I realise the broadness of the questions to further wittle down my queries from the already useful answers,

How do I create such a cloud library of modis data ?
How is such data then further automated using Python, specifically are there any examples of such automation or any places I need to go to learn this specific parts of Python?

If this is still too basic I apologise, but I am coming at this from a complete Python beginners perspective.