texts: Calculate and plot the feature importance of the input dataset of a random forest classification within Google Earth Engine

texts:
I run a random forest classification for agricultural land use and other land cover classes (12 classes). My dataset for 2019 consists of a Sentinel 1 and Sentinel 2 monthly time series, statistical phenometrics (e.g. min/ max/ median/ stdDev) and ancillary data as texture metrics and topographic derivates. I included a lot of different indices calculated on the basis of the sentinel 2 optic bands. This resulted in a single image with 294 bands as a big input data cube for the random forest algorithm. In order to decrease computational time I would like to calculate the feature importance of the different bands in that image for each one class against all classification and then only use the most important bands. So far I didn't find a lot in different blogs and platforms.
I found an approach on MyGeoBlog that is unfortunately showing an error that I couldn't resolve.  In theory this approach should work with the one class against all other classification approach. Find my adjusted script with the following link (note that its a different dataset with only 216 bands without the sentinel 1 data):
https://code.earthengine.google.com/?scriptPath=users%2Fmarcellitomohr93%2FArmenia_LULCC%3ALULC_Classification2019%2FGoris2019_ClassificationIterations%2FFeatureImportance_GrasslandMask
The original script on MyGeoBlog you may access via this link (however it is lacking input data):
https://mygeoblog.com/2020/02/12/random-forest-trees-and-variable-importance/
In my case GEE shows the usual point feature size exceeding error but also the error: "Error generating chart: Property 'Grass_Mask' of feature '0' is missing." Since for this script you can only calculate probability with 2 classes I transformed my 12 classes of a LULC_Class property to a binary property column that is called grass mask (with only values of 1 for grassland and 2 for all other classes).
I'm happy to hear about your thoughts and experiences with calculating the feature importance in the GEE environment. I could also do it in R but then I have to download my very large dataset.


texts: As I loaded sDNA as a plug-in in ArcGIS 10.6 following the instructions, I found most of the tools invalid as the screencut shows. And the tools that seem loaded normally, like the network preparation report error when I load a network file. Are there any problems with my installation?


texts: Location allocation Road Block- Unable to "Solve"

texts:
Software: ArcGIS/ArcMap 10.3.1 with network analyst extension
Objective: Perform location-allocation services, based off of two Required Facilities and Demand Points.

I feel like I have all critical pieces I need:

State Shapefile, created a network dataset via it (most recent from here: http://download.geofabrik.de/north-america/us/florida.html)
Facilities Locations - Created a Shapefile (Through ArcCatalog .xls to shapefile) Properties changed to make them Required Facilities, Respective Facility capacities entered
Demand Points - Created a shapefile (through ArcCatalog .xls to shapefile)
Location Allocation- Accepts my network dataset, impedance cut-off at 120, Facilities to choose at 2, Problem type changed to Maximize Coverage, and Travel Point selected "Demand to Facility".

However, I believe I have a layer connectivity issue. I say this because when I attempt to "Solve" my analysis. No lines are rendered. I've tried countless different file types/sources (shape, osm, pbf, bz2), as well as trying to allow ArcMap Geocode my addresses, importing my LAT/LONG with their X/Y fields, with no avail.

Everything on my map, looks great in terms of the roads, where the facilities and demand points get plotted. However, when attempting to "Solve" my analysis... no lines are drawn.

Does anybody have suggestions of what I could possibly be doing wrong?


texts: Tools to export shps from Oracle database

texts: What open source tools could be used to export shapefiles from an Oracle database. I have tried using QGIS but it does crash way too much.


texts:
I am using the GDAL Java bindings to load vector data sets into a Java based GIS application.
The GDAL Java bindings seem to contain no classes or methods to parse the Feature#getStyleString into something I can easily work with. Looking at the native GDAL code, there exist classes like OGRStyleTool and enumerations like ogr_style_tool_class_id which do the parsing for you.
What are my Java options if I want to extract the styling information from the style string ? Do I have to write my own parser for this, or am I overlooking another option ?


texts:
I am trying to show some elevation/cross section data in a graph. I am trying to do this in arcgis graphs. But because of my previous graph work was in .net, I found arcobjects.net API for graphing is not so much coder friendly. So, now I want to go back in .net graphs. But before that I want to be confirm that, I am not doing any kind of mistakes. So I goggled for advantage and disadvantage of arcgis graphs from coder perspective. But I couldn't find that much information there. So I want some help from experts.

EDIT :
Based on selecting water-body from the map, I am showing a graph. This graph shows elevation value of that water-body's elevation. If user selects another water-body, then another series will be added. I am showing this graph from my add-in. It's not necessary to show it within ArcGIS Desktop's graphing framework. But if I can show it from ArcGIS Desktop's graphing framework, it will be great.

For graphing, which one would be good .net graphs or arcgis graphs?


texts: PyQGIS is working too slow. Can I use GeoPandas?

texts:
I have a layer with classification od buildings. I want to add "Area" to my attributes table and I use PyQGIS. For another layer with 50 records it works, but my layer is huge, contains 6.8mln record and it is working too slow. With another problem I used GeoPandas which is working so fast. I will paste my code and I need advice. Can I do it with GeoPandas?
from qgis.core import *
from qgis.utils import *
from qgis.analysis import QgsNativeAlgorithms
from PyQt5.QtCore import QVariant
from qgis.core import QgsApplication, QgsProcessingFeedback, QgsRasterLayer
import sys
import geopandas

sys.path.append('/usr/lib/qgis')
sys.path.append('/usr/share/qgis/python/plugins')
os.environ["QT_QPA_PLATFORM"] = "offscreen"

QgsApplication.setPrefixPath(r'/usr/bin/qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()

import processing
from processing.core.Processing import Processing

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
feedback = QgsProcessingFeedback()

#adding the Area field
layer = QgsVectorLayer(r'/home/gis/polskagisencoding.shp', "polskagisencoding", "ogr")
provider = layer.dataProvider()
area_field = QgsField("Area", QVariant.Int)
provider.addAttributes([area_field])
layer.updateFields()

#updating the Area field for each feature
idx = provider.fieldNameIndex('Area')
for feature in layer.getFeatures():
    attrs = {idx : int(feature.geometry().area())}
    layer.dataProvider().changeAttributeValues({feature.id() : attrs})


texts: Subtract two pixels-only- layers using rectangle extent

texts:
I have coe in GEE that calcultes the NDVI in a given polygon, from 2000 (based on landsat 7) and 2019 (landsat 8) ,create for each one a layer with pixels that their NDVI value is higher than 0.3 and then subtract them from each other in order to detect change.

The problem is that when I run my code for some reason the result is that ther is no diffference and I believe I have problem in my code.
I think mabybe GEE doesn't do the subtruction if the pixels are not on the same location so I always get that no change.

maybe it will be easier to understand from my code:
https://code.earthengine.google.com/98102459e86cc67232896d226f19565b

Basically I have create to layers in a long process and in this specific line:

// Map.centerObject(geometry,10);
for (var i in listOfNumbers1){
  var imageNDVI1 = ee.Image(listOfImages.get(listOfNumbers[i]));
  var hot2019=imageNDVI1.gt(0.3)
  .selfMask()
  .rename('hotspot03');
Map.addLayer(hot2019,{palette: '0000ff'},'2019');}

and then I subtract the two layers (hot2000 and hot 2010) like this:

var difference=hot2019.subtract(hot2000);

Map.addLayer(difference,{min:-1, max:1 ,palette: ['ff3300', 'a3a375','009933']},'difference');

My end goal is to get new layer that will be calssified into 3 values: -1,0 and 1

*I thought maybe the reason it happened is because it doesn't have rectangle extent
*Any polygon would work. I have checked Delhi and New Mexico


texts: QGIS raster image transparency, NO DATA value for all raster images?

texts:
I'm working with sets of aerial images in QGIS which are geocoded including orientation i.e. rotated.

After loading images QGIS draws a black box around each image that extends from the far west to far east pixel and from the far north to far south pixel. Within this black box lays the actual image in it's correct orientation.

To get rid of this black box I currently enter 0 in Layer Properties / Transparency / Additional NODATA value but the change has to be made on every single image. Unfortunately the transparency setting is not part of the Style definition and therefore cannot be copied from one to all other raster layers / images.

Is there better way to add a global NODATA value to all raster data in current workspace?

Raster Calculator won't work as there are no NODATA areas in the actual image files.


texts: OpenLayers does not display vector tiles in their origin projection

texts:
I have 3 vector layers and one raster which I want to display with OpenLayers.
I used Geoserver and Geowebcache for creating VectorTiles cache (in EPSG:4326). These cached  layers displayed well in the standard Geoserver viewer.

Then I created map in OL (for this example I left only one raster (basemap) and one vector tile layer):

var layerCity = 'xxx:city',
    layerTerrain = 'xxx:terrain';

var epsg_projection = '4326';

var City = new ol.layer.VectorTile({
    source: new ol.source.VectorTile({
      tilePixelRatio: 1, // oversampling when > 1
      tileGrid: ol.tilegrid.createXYZ({maxZoom: 13}),
      format: new ol.format.MVT(),
      projection: 'EPSG:4326',
      url: 'http://xxx/geoserver/gwc/service/tms/1.0.0/' + layerCity + '@EPSG%3A' + epsg_projection +'@pbf/{z}/{x}/{-y}.pbf'
    }),
    style: function(feature) { //sdsds
      return styleCity = new ol.style.Style({
              image: new ol.style.RegularShape({
                  stroke: new ol.style.Stroke({
                    color: 'white',
                    width: 1.5
                  }),
                  fill: new ol.style.Fill({
                    color: 'black'
                  }),
                  radius: 5,
                  points: 4,
                  angle: Math.PI / 4,
             }),
             text: new ol.style.Text({
                stroke: new ol.style.Stroke({
                   color: '#ffffff',
                   width: 2.5
                 }),
                // fill: new ol.style.Fill({
                //    color: 'f2cb00',
                //    width: 2.5
                //  }),
                font: 'bold 14px Roboto, sans-serif',
                text: feature.get("name").toString(),
                offsetX: 17,
                offsetY: -12
            })
        })
    }
  });

var Terrain = new ol.layer.Tile({
    source: new ol.source.TileWMS({
      url: 'http://xxx/geoserver/gwc/service/wms?',
      params: {
        'LAYERS': layerTerrain,
        'TILED': true,
        'VERSION': '1.1.1',
        // 'SRS': 'EPSG:4326',
        'FORMAT': 'image/png8'
      }
  })
});

var map = new ol.Map({
    target: 'map-frame',
    view: new ol.View({
        projection: 'EPSG:4326',
        // center: [-26.103, -4.195], //Deois
        zoom: 2,
        minZoom: 0,
        maxZoom: 8
      //  projection: 'EPSG:4326'
    }),
    layers: [Terrain, City]
});

As you сan see I set projection for my map to "EPSG:4326", because default prj is 3857.

But unfortunatly I got a partial result: displaying of basemap without vectortile layer visualisation. Also I used Chrome developer extension for checking donwloads. So I got about 15 pbf files without fails but I don't see that on the map!

However, if I change (or remove) map projection to 3857, I get the vector tile visualisation. I think it must sounds strange because on the server I have tiles' directories only with such format: "EPSG_4326...".

I have no any idea for solving this problem.


texts: Determine which polygons are cut by circumferences

texts:
In QGIS, I have a set of points from which I created a series of voronoi polygons. I also have a set of circumferences created from another points, like this:

I need to identify which voronoi polygons are cut by which circumferences (one polygon would appear in more than one circunference, so I need to track de attributes), but everything I have seen is about intersection.

Bonus points would be polygons with only area > N, but I need to start somewhere.

EDIT:

First solution Joseph gave me in the comments gave me the polygons which intersected with the circumferences as a whole, in a similar fashion to wiltomap. The actual answer from Joseph worked for me.

As an extra, from the anwer I was able to calculate the polygon area which allows me to filter those with a very small one.


texts: Can legends be exported as image in CARTO BUILDER?

texts: How to add a legend when exporting a BUILDER map as jpg/png image?


texts: Add SLDs to GeoServer using REST API Fails

texts:
I have been writing some python code to automatically add data layers from PostGIS to GeoServer and attach a suitable SLD as default. I got this working using the following code:

store=("<dataStore><name>data</name><connectionParameters><host>localhost</host><port>5432</port><database>case</database><schema>public</schema><user>" + postgisusername + "</user><passwd>" + postgispassword + "</passwd><dbtype>postgis</dbtype></connectionParameters></dataStore>")

c = pycurl.Curl()
c.setopt(c.VERBOSE, 1)
c.setopt(c.HTTPHEADER, ['Content-type: text/xml'])
c.setopt(c.USERPWD, username + ":" + password)
c.setopt(c.URL, 'http://localhost:8080/geoserver/rest/workspaces/test/datastores')
c.setopt(c.CUSTOMREQUEST, "POST")
c.setopt(c.POSTFIELDS, store)
c.perform()
c.close()

sldname='<style><name>DataLine</name><filename>DataLine.sld</filename></style>'

c = pycurl.Curl()
c.setopt(c.VERBOSE, 1)
c.setopt(c.HTTPHEADER, ['Content-type: text/xml'])
c.setopt(c.USERPWD, username + ":" + password)
c.setopt(c.URL, 'http://localhost:8080/geoserver/rest/styles')
c.setopt(c.CUSTOMREQUEST, "POST")
c.setopt(c.POSTFIELDS, sldname)
c.perform()
c.close()

filesize = os.path.getsize("Data Line.sld")
f = open("Data Line.sld", 'rb')

c = pycurl.Curl()
c.setopt(c.VERBOSE, 1)
c.setopt(c.HTTPHEADER, ['Content-type: application/vnd.ogc.sld+xml'])
c.setopt(c.USERPWD, username + ":" + password)
c.setopt(c.INFILESIZE, filesize)
c.setopt(c.URL, 'http://localhost:8080/geoserver/rest/styles/DataLine')
c.setopt(c.PUT, 1)
c.setopt(c.INFILE, f)
c.perform()
c.close()
f.close()

sldtolayer='<layer><defaultStyle><name>DataLine</name></defaultStyle><enabled>true</enabled></layer>'

c = pycurl.Curl()
c.setopt(c.VERBOSE, 1)
c.setopt(c.HTTPHEADER, ['Content-type: text/xml'])
c.setopt(c.USERPWD, username + ":" + password)
c.setopt(c.URL, 'http://localhost:8080/geoserver/rest/layers/dataline')
c.setopt(c.CUSTOMREQUEST, "PUT")
c.setopt(c.POSTFIELDS, sldtolayer)
c.perform()
c.close()

This works very well for all SLDs where I am not referencing an external graphic such as an SVG for symbology. I moved a folder of SVGs into the styles folder within the data directory and restarted the OpenGeo Suite, to make sure it knew the SVGs existed.

However, the code creates a blank SLD within GeoServer.

The error message I receive is a Java Unknown Source error and external graphic location.

I know this SLD works if I load through the GUI, so it must be something going on through the REST API.

Can anyone help?

thanks


texts:
I am using QGIS v1.6 on Ubuntu 10.10.

When I use the Fetch Python Plugins menu item QGIS crashes. If I run from the terminal with sudo I get the following error:

Segmentation Fault

Does anyone know how to fix this? I have this problem on 2 installs of Ubuntu on 2 different computers.


texts: Averaging a big number of pixels within a several km radius

texts: I work with ArcMap software and now I have the aim to average individual pixels within a several km radius from daily data for at least 1 month. So I made the grid with 4 km step (the step size is not important for the question) and now I have to average pixels according to the described method. How do you think what the tool from ArcTools window I should use for this?


texts: Need to develop a model for measuring inter-annual and seasonal variability in pasture cover using modis 16 day composites

texts:
Aiming to use statistical methods in ArcGis and Envi image analysis to measure inter-annual and seasonal variability in paster cover using MODIS 16 day composite image from 2000-2014.

Quite stuck on where to start and need any information or links to papers or previous analysis' of a similar style.


texts: GGeoXml error with Arcgis extension for google maps api

texts:
I'm using the ArcGis google maps extension.

I've created a simple page based on the esri example here...

I have successfully buffered a polyline, however, when I create a query and execute the related queryTask, I see the following error...(in Firebug)

"GGeoXml object could not load the KML file correctly"

as a rough guess, I'd say there's a problem with my proxy setup, anyone have any ideas?


texts: ArcGIS Pro missing records in attribute table from SQL Server?

texts:
I created a table to test ArcGIS Pro 2.2.3 with an Azure SQL Server. The structure of the table is simply an OBJECTID, a shape (Point) and a text field. The table contains a set of 2 points with different locations.

I am then loading this table as a Feature Class into a map in ArcGIS Pro. I can see the 2 points and I can select them both from the map view. The issue I have is that the Attribute Table only lists a single point when listing all records. Furthermore I get the following behaviour from the Attribute Table:

If looking at all records with nothing selected on map: I see 1 record and caption says "0 of 1 selected"
If looking at selected records with nothing selected on map: I see no record and caption says "0 of 1 selected"
If looking at all records with both selected on map: I see 1 record and caption says "2 of 1 selected" (see image below)
If looking at selected records with both selected on map: I see 2 records and caption says "2 of 1 selected"

The table itself was not created by ArcGIS Pro but was defined in SQL as a table with a spatial field using the following code. The database is hosted on an Azure SQL Server:

/****** Object:  Table [geotest].[geometry] ******/
SET ANSI_NULLS ON
GO

SET QUOTED_IDENTIFIER ON
GO

CREATE TABLE [geotest].[geometry](
    [OBJECTID] [int] IDENTITY(1,1) NOT NULL,
    [Shape] [geometry] NULL,
    [col_as_text] [nchar](100) NULL,
CONSTRAINT [PK_geometry] PRIMARY KEY CLUSTERED
(
    [OBJECTID] ASC
)WITH (PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON)
)

GO

/****** Object:  Index [FDO_Shape] ******/
CREATE SPATIAL INDEX [FDO_Shape] ON [geotest].[geometry]
(
    [Shape]
)USING  GEOMETRY_AUTO_GRID
WITH (BOUNDING_BOX =(-400, -90, 400, 90),
CELLS_PER_OBJECT = 16, PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, SORT_IN_TEMPDB = OFF, DROP_EXISTING = OFF, ONLINE = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON)
GO


texts: How to open and display "complex type" GeoTIFF?

texts:
Here is a GeoTIFF (single-band) that could be successfully loaded into QGIS but produces no visual output regardless its style settings/parameters. The file can be found here here. Other GIS (e.g., ENVI) could load and display the same file with no problem. What is the problem with this GeoTIFF file? Is it in a format that is currently not supported by QGIS?
[Original question]
The "complex type" GeoTIFFs cannot be opened/displayed by QGIS.
The image is opened without error message but nothing is displayed. In addition, all the pixels values of the image are indicated as "nan" by the "Identify feature" tool.
This GeoTIFF type is characterised by the fact that every pixel is a complex value.
This type is handled properly by Python/GDAL, IDL/ENVI...
Similarly to ENVI, the absolute value of every pixel should be displayed, and the "Identify feature" tool should indicated the couple of Real and Imaginary values related to the pixel clicked.
Please find here an example of complex GeoTIFF that can be properly displayed with other GIS softwares (like ENVI).
Handling this file type with QGIS would be very useful as many SAR images are stored with complex values.


texts:
I am attempting to build an ArcGIS Python script within a basic Toolbox so I can scan through and clip all shapefiles within predefined folders.

This is easy with the Geoprocessing tool "Clip" for single features or multiple with using batch. However, I have reached a very frustrating impasse and haven't been able to find a solution.  I wrote a basic code below to show the problem I am having... what am I doing wrong?

import arcpy, os
arcpy.AddMessage("Basic Clip Attempt")
Clip_Feautre = arcpy.GetParameterAsText(0)
# "X:/3 GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Boundaries/Boundary_Extended"
Target_Feature = arcpy.GetParameterAsText(1)
# "X:/3-GIS/5_BaseData/AltaLIS/AltaLIS_Hydrography_160619/BF_HYDRO_POLYGON.shp"
Save_Location = arcpy.GetParameterAsText(2)
# "X:/3-GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Scratch/Hydro_Scratch"

if os.path.exists((Target_Feature)):
    arcpy.AddMessage(Target_Feature + "----- Path Exists")
else:
    arcpy.AddMessage(Target_Feature + "----- Path Does Not Exist")

# Guess what... the file exists!
# X:\3-GIS\5_BaseData\AltaLIS\AltaLIS_Hydrography_160619\BF_HYDRO_POLYGON.shp----- Path Exists
try:
# The following is a copy from an ArcGIS geoprocessing run that worked.
arcpy.Clip_analysis(in_features="X:/3-GIS/5_BaseData/AltaLIS/AltaLIS_Hydrography_160619/BF_HYDRO_POLYGON.shp", clip_features="X:/3-GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Boundaries/Boundary_Extended", out_feature_class="X:/3-GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Scratch/Hydro_Scratch", cluster_tolerance="")
except:
arcpy.AddMessage(arcpy.GetMessages(2))
# ERROR 000229: Cannot open X:/3-GIS/5_BaseData/AltaLIS/AltaLIS_Hydrography_160619/BF_HYDRO_POLYGON.shp
# Failed to execute (Clip).

arcpy.AddMessage("Attempt 1 has failed.  Attempt 2 in progress")

in_features_raw = r'X:/3-GIS/5_BaseData/AltaLIS/AltaLIS_Hydrography_160619/BF_HYDRO_POLYGON.shp'
clip_features_raw = r'X:/3-GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Boundaries/Boundary_Extended'
out_feature_class_raw = r'X:/3-GIS/1_Projects/2017_Bacardi_Creek/06_Geodatabase/2017_Bacardi_Creek.gdb/Scratch/Hydro_Scratch'
try:
    arcpy.Clip_analysis(in_features_raw,clip_features_raw,out_feature_class_raw,cluster_tolerance="")
except:
    arcpy.AddMessage(arcpy.GetMessages(2))
    # ERROR 000732: Input Features: Dataset X:/3-GIS/5_BaseData/AltaLIS/AltaLIS_Hydrography_160619/BF_HYDRO_POLYGON.shp does not exist or is not supported
    # Failed to execute (Clip).
    arcpy.AddMessage("Attempt 2 has failed... Attempt 3 in progress...")
    try:
        arcpy.Clip_analysis(Target_Feature,Clip_Feautre,Save_Location,cluster_tolerance="")
    except:
        arcpy.AddMessage(arcpy.GetMessages(2))
        # ERROR 000229: Cannot open X:\3-GIS\5_BaseData\AltaLIS\AltaLIS_Hydrography_160619\BF_HYDRO_POLYGON.shp
        # Failed to execute (Clip).
        arcpy.AddMessage("Attempt 3 has failed... Attempt 4 in progress........")
        try:
            Clip_Feautre_replace = Clip_Feautre.replace("\\","\\\\")
            Target_Feature_replace = Target_Feature.replace("\\","\\\\")
            Save_Location_replace = Save_Location.replace("\\","\\\\")
            arcpy.Clip_analysis(Clip_Feautre_replace,Target_Feature_replace,Save_Location_replace,cluster_tolerance="")
        except:
            arcpy.AddMessage(arcpy.GetMessages(2))
            # ERROR 000229: Cannot open X:\\3-GIS\\5_BaseData\\AltaLIS\\AltaLIS_Hydrography_160619\\BF_HYDRO_POLYGON.shp
            # Failed to execute (Clip).
            arcpy.AddMessage("Attempt 4 has failed... .... I don't know...")


texts: Geoserver 2.4.5 with Openlayers 2.13.1 overlay WMS layer

texts:
I'm using GeoServer with OpenLayers, happens that my WMS layers do not appear.
Do I need to overlay a WMS layer in PostGIS that I have with google maps api and does not work appears in WMS layer switcher but does not appear on the map. Generated the wms link to the Earth and visualized layer is perfect.
I then tested by creating a simple map in openlayers with those with a WMS layer that comes in geoserver and neither worked.
What could be wrong? Debuguei the code in chrome, with a breakpoint in loading the WMS layer and does not give any errors.

This is my code:

<!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>

 <script src="http://maps.google.com/maps/api/js?sensor=false&v=3.2"></script>

 <script type='text/javascript'>

 var map;
 var mercatorProjection = new OpenLayers.Projection('EPSG:900913');
 var latLongProjection = new OpenLayers.Projection('EPSG:4326');
 var bounds = new OpenLayers.Bounds(
                -51.2809219360352, -30.2445888519287,
                -51.0207977294922, -29.9661273956299
            );

 function init() {

  map = new OpenLayers.Map('map_element',{
    maxExtent: bounds,
    maxResolution:  0.0010877400636672,
    units: 'm',
    //allOverlays: true,
    projection: mercatorProjection,
    displayProjection: latLongProjection
 });

 var google_streets = new OpenLayers.Layer.Google(
    "Ruas",
    {numZoomLevels: 20}
    );

var wms_layer = new OpenLayers.Layer.WMS(
'Linhas Ônibus',
'http://localhost:8080/geoserver/Teste/wms',
{layers: 'Teste:onibus', transparent: true},
{isBaseLayer: false,
    opacity: 0.7}
);

//Adiciona as camadas ao mapa
map.addLayers([google_streets, wms_layer]);

var point = new OpenLayers.LonLat(-51.22,-30.08);
point.transform(new OpenLayers.Projection("EPSG:4326"),
    map.getProjectionObject());
map.setCenter(point, 10);
 //Camada de controle que vai mostrar as camadas no mapa
 map.addControl(new OpenLayers.Control.LayerSwitcher({}));

 //Verifica se o mapa tem um ponto central e o extende a sua extensão máxima
 if(!map.getCenter()){
 map.zoomToMaxExtent();

 }
 }

 </script>
</head>

<body onload='init();'> <!-- Chama a função js init() -->

<!--Elemento HTML onde o mapa é exibido -->
 <div id='map_element' style='width: 800px; height: 800px;'>
 </div>

</body>
</html>


texts: Updating geometry as Multipolygon from Polygon table using PostGIS?

texts:
There are two tables in PostGIS. In which, one has polygon geometries with layer name wise and another one table has master layer name list.

The requirement is to collect the polygon geometries as layer wise from first table and update it as Multipolygon in the geometry column of second table that is layer master table. Below screenshot describes the requirement.

But boundaries of the each polygon should be displayed when viewing the second table on QGIS/Pgadmin graphical viewer after updating them as Multipolygon.

Seeking a sample PostGIS queries to achieve the above.


texts:
I'm trying to obtain the same label style than Perpendicular Placement Orientation from ArcGIS in QGIS. As QGIS doesn't have this feature as default, the option would be using custom placement for this.
An ArcGIS example:

A workaround in QGIS
Settings:
Data defined coordinate X (label in the middle of he line):
$x_at( num_points( $geometry)/2)

Data defined coordinate Y (same as X's location):
$y_at( num_points( $geometry)/2)

Data defined Rotation:
angle_at_vertex( $geometry , num_points( $geometry)/2)

For only for 2/3 labels, the placement is acceptable. For lines with a lot of breaks with different angles, it's more difficult to obtain a nice label.
ArcGIS seems to search the line segment with fewer breaks or a straight segment. Is there a way to obtain coordinates or geometry of the straightest line segment using Expression builder for a better label placement?


texts: Routing for specific subjects (cars, pedestrians, ..) dynamically

texts:
I need create app which will be provide routing for specific category. For example only for cars, only for bikes or only for pedestrians. For importing OSM data I'll used default mapconfig for all way types.

Can you tell me which classes are officialy used for cars, pedestrians and how?

I wants use it into pgr_dijsktra function where I will be able select specific ways.

I'll tried this but not works. When I tried same script without where condition for classes it works, but I need specific category dynamically in condition.:

SELECT e.seq, e.path_seq, e.node, e.cost, e.agg_cost, w.name, w.the_geom, w.x1, w.y1, w.x2, w.y2 FROM pgr_dijkstra('
SELECT gid AS id,
     source,
     target,
     cost_s AS cost,
     reverse_cost_s AS reverse_cost
    FROM pgr.ways where class_id in (101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116)',
    35959,
    19303,
directed := true)e
left join pgr.ways w on w.gid = e.edge
order by e.seq

How can I do that?


texts: How to model money cost as an accumulator in network analysis?

texts:
This is not necessarily a technical coding question, although any links that include coding are welcome if they are answering my question.

I am looking for research or analysis that has been done in network analysis using the money cost of travel as a cost accumulator per segment of road/rail/bus route/etc. Does anybody know of a good resource such as articles/publications/blogs/journal articles showing someone who either does it, suggests how to do it, or explains the methodology behind it?

I am looking into the possibility of using money cost as a cost accumulator in the ArcGIS 10.0/10.1 network analyst least-cost path problem. My thoughts are, if a commuter rail were to be built between two metropolitan centers, which work commuters would be likely to choose the commuter rail over driving a car the whole way when they are factoring in how much it costs them to make the trip?


texts: Using variable from arcpy.GetParameterAsText() as part of output path?

texts:
I´m trying to use a variable as part of the output path of a shutil command.
Unfortunately I couldn´t work it out how to do it. Can anybody help?

import arcpy
import shutil
Name = arcpy.GetParameterAsText(0)
shutil.copyfile("C:\\test.txt", "C:\\%Name%.txt")

This question here is similar, but was never answered:
Using ArcGIS ModelBuilder to perform In-line variable substitution for input data path?


texts: Transforming a GRASS formula into a code usable in R

texts:
I would like to transform the formula at the botton into a formula for R , as R runs better and fastly with files I have, but I am new with it- For me it is enough even a "literary" traduction, where I can start from, or even some tips. I try to explain the formula below, it is a calculation of an extreme climate event projection index, the Warm Spell Duration Index (WSDI), inside species range maps (around 5000 mammals).
DEFINITION. It is a "modified" WSDI, calculated as follow: annual count of days, WITHIN EVERY SPECIES RANGE, with at least 6 consecutive days, centered on a five days window, when Tmax > 90th percentile of the historical ("1961-1990 average period" in my case). So, differently from a global index, the 90th percentile is different for each species and calculated, spatially, once again, for every iteration of external loop (as saying, for every species range map).
FORMULA EXPLANATION. It is a nested loop of: WSDI calculation (internal loop) for all mammals species (external loop).I use 'g.region' and 'r.mask' to work within species range, and then use the 'eval' command to extract 90th percentile from the unvariate statistic of historical raster. WSDI is calculated for year 2091, so from raster #1827 to raster #2191 in the data folder I have of rasters of daily Tmax projections (the internal loop). Line 10 (which is a bit long), with the help of temporary raster ('temp_WSDI_rcp26_2091') and 'max' function (line 11) is the basis for WSDI calculation (6 consecutive days when Tmax > 90th percentile of the historical). 'WSDI_"$species"_rcp26_2091' at the end of every internal loop is the final WSDI*Range for every species. Finally I remove mask and end external loop.
DIFFICULTIES. Main difficulties are in creating raster loop (the external loop), as it uses raster instead of numbers (as the internal loop use numbers), and also in setting 'region' and 'mask' (both allow in confining 90th percentile calculation within single species range). 'calc' should be used instead of 'r.mapcalc', but I still was not able to understand how to define different 'function' ( from R manual: calc(x, function, filename='', na.rm, forcefun=FALSE, forceapply=FALSE, ...) ). The 'if' command should be the same, maybe changing language. Finally, I think 'eval' is a bash command, so kind of "universal" for all software, hoping I can use it in this case. Here is the formula:
1. for species in `g.list rast pattern=* mapset=Range_M`;
2. do
3. g.region raster="$species"@Range_M;
4. r.mask raster="$species"@Range_M;
5. eval `r.univar -ge  map=1961_1990_tasmax_historical`;
6. r.mapcalc expression="temp_WSDI_rcp26_2091=0" --overwrite;
7. r.mapcalc expression="WSDI_"$species"_rcp26_2091=0";
8. for i in {1827..2191};
9. do
10. r.mapcalc expression="temp_WSDI_rcp26_2091=if(tasmax_day_MRI_CGCM3_rcp26_r1i1p1_20860101_20951231.$i-2>$percentile_90 && tasmax_day_MRI_CGCM3_rcp26_r1i1p1_20860101_20951231.$i-1>$percentile_90 && tasmax_day_MRI_CGCM3_rcp26_r1i1p1_20860101_20951231.$i>$percentile_90 && tasmax_day_MRI_CGCM3_rcp26_r1i1p1_20860101_20951231.$i+1>$percentile_90 && tasmax_day_MRI_CGCM3_rcp26_r1i1p1_20860101_20951231.$i+2>$percentile_90,temp_WSDI_rcp26_2091+1,0)" --overwrite;
11. r.mapcalc expression="WSDI_"$species"_rcp26_2091=max(WSDI_"$species"_rcp26_2091, temp_WSDI_rcp26_2091)" --overwrite;
12. done;
13. r.mask -r;
14. done


texts: Joining a csv table to the attribute table of a feature class in arc gis 10.3

texts:
I am new to arcpy. I was trying to join a csv table to the attribute table of a feature class. Here I would like to share what i was trying to do.

import arcpy, csv

from arcpy import env

env.workspace = r"F:\Users\User\Hydrology\Test.gdb"

fcl = arcpy.ListFeatureClasses("*")
for fc in fcl:

 print "Joining table to " + fc

HRU06 = "F:\Users\User\River\HRU06.csv"

HRU08 = "F:\Users\User\River\HRU08.csv"

HRU10 = "F:\Users\User\River\HRU10.csv"

HRU12 = "F:\Users\User\River\HRU12.csv"

arcpy.JoinField_management("catchments_03", "FID", HRU06, "HRU_06_FID")
arcpy.JoinField_management("catchments_04", "FID", HRUC08, "HRU_08_FID")
arcpy.JoinField_management("catchments_05", "FID", HRU10, "HRU_10_FID")
arcpy.JoinField_management("catchments_06", "FID", HRU12, "HRU_12_FID")

But I am not getting it done. Can someone please help me? I want to join the csv files on the basis of "FID" field of the catchments.


texts: Unselect all layers and clear the legend in Openlayers

texts:
Is the there a way of clearing the legend in Openlayers by unselecting all the layers? We have a map that starts with several layers.

I was thinking that we might be able to add a button to the map that would do this.


texts: Using Python parser of ArcGIS Field Calculator to populate values in one field based on those in another?

texts:
I am trying to populate a field with a sorting code that represents the value of another field. Example:

 EXPCODE|SORT|
 S4DNN2 |3.5
 H2BNN2  5.0

EXPCODE is a field with a value like "S4DNN2." SORT is a field that needs to represent the EXPCODE with a numerical value. What would I use for Python code to automatically populate the SORT field? I started with a simple "if" statement, but can't quite figure it out.

def EXPCODE(field_value):

    if field_value = "S4DNN2"
        return field_value "3.5"

    if field_value = "H2BNN2"
        return field_value "5.0"

The EXPCODE values I show are examples, I have about 67 different ones with unique numerical values that need to be populated in "SORT"

Any ideas?

Update: I am trying to execute this with the field calculator with use of the python parser


texts: How to avoid multiple "Other" Columns with similar Styling in Carto

texts:
I have a shapefile for all US ZCTAs. Each ZCTA has been clustered by proximity and given a value in a new column called Zip_Rollup.

I am trying to add a style by value using the Zip_Rollup values.  This works great except that it only allows 11 unique values and the rest are lumped as “others” with the same color making it hard to distinguish them as unique regions.

Zip_Rollup has 1000 + unique values which I would like to visualize.

It’s fine if there are duplicate colors. Is there a way to apply a random color to polygon fill by value in Carto, or other option to solve this problem?

I found a similar question with answer for Cartodb at the link below, but could not translate this to Carto myself.

Automatically coloring a large number of polygons in CartoDB


texts:
I'm updating a cached map service by features at different scale levels to only recreate tiles over changed data.

It works fine, but it occurred to me it might be more efficient to cache by super-tile extents, so if multiple of my caching features are within the same super-tile area it doesn't re-cache the same area multiple times.

Does anyone know if there is a way to visualize or map super-tile areas?

I know I can see the bundle areas with the task status fc in the cache status.gdb but the bundle areas are potentially much larger than a super-tile.

I'm using ArcGIS Server 10.1 sp1.


texts:
I'm looking for Image datas for East Africa. Specifically for Rwanda. Does anyone know if such data are available at all free? My target is to digitize roads and buildings with the help of image datas.

I am working with Grass GIS. My job is to find small hydro Power sites. I already did that. For these locations special boundary conditions must be accomplish. E.g. Accessibility of a location (maximum distance to roads about 2000 meter) and a maximum distance to the consumer needs to be defined. Because there are no vector datas of buildings and roads, I need the Image Datas to digitize them.


texts:
I need to use the Lambert Conformal Conic for Natural Resources Canada EPSG:42304, which I could not find in QGIS 2.10.1 CRS list. I located the proj4 parameters on-line from three locations, all of which were the same:

+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs

But when I add this and test it in the Custom Coordinate Reference System Definition a pop up comes with a warning saying 'Northing and Easting must be in decimal form'.

I have tried looking around for an answer or a different set of parameters, but to no avail.

Can anybody help me?


texts: Displaying a device which has changed it's connection ID in the last week using Arcade Expression in ArcGIS Online?

texts:
I have about 250 devices in a map, but I want to only display devices that have had a certain attribute change within the past week.

I have minimal experience with Arcade expressions and cannot find guidance for this particular query online.


texts: Downloading Subarea vs Clipping Images to AOI using Planet API?

texts:
This question is associated with the two Planet API tutorials: "Download a Subarea" and "Clip Images to AOI".

It seems to me that both tutorials demonstrate methods to download only a piece of an image from a large scene so as to increase efficiency.

How do these tutorials differ?


texts: Simulate a click event on FeatureLayer in ESRI JavaScript API

texts:
When my app changes context (user clicks on a tab) I want to programmatically "Select" the first of 6 points in my Feature Layer so I get the blue box around the point and activate the info window.

I looked around the API Docs and found the FeatureLayer's selectFeatures method.  This method takes a Query object, and a selection method (I'm creating a new selection so I'll use SELECTION_NEW).  The function returns a promise, but do I really need any kind of callback? I'm just trying to select a feature.

Here is the definition of my feature layer:

var projectsAGOL = new esri.layers.FeatureLayer("http://services.arcgis.com/QVENGdaPbd4LUkLV/arcgis/rest/services/rioLocoSites/FeatureServer/0",{
    infoTemplate: infoTemplate,
    visible: true,
    outFields: ["*"]
});

Click handler for the Tab:

jQuery("#prjMap").click(function(){
    var query = new esri.tasks.Query(); // Create a Query Object
    query.where = "OBJECTID = 0"; // Select the first point from the Feature Layer
    projectsAGOL.selectFeatures(query, esri.layers.FeatureLayer.SELECTION_NEW); // Apply selection
});

This code does not produce the desired outcome nor does it produce an error.  Nothing happens.  Am I even using the correct method to make a selection?  I'm more familiar with the Leafletjs mapping library.

Update: Provided jsFiddle


texts:
If I use a single WMS layer as map base layer it loads the map, but adding different layers like (Google, Yahoo) to map it won't work. Can anyone tell me how to overcome this...
please check this:

    <script src="http://maps.google.com/maps/api/js?sensor=false&v=3.3" type="text/javascript"></script>
    <script src="http://api.maps.yahoo.com/ajaxymap?v=3.0&appid=abcdefg1234" type="text/javascript"></script>
    <script src="http://openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
    <script type="text/javascript">
        var map;
        function init() {
           map = new OpenLayers.Map('map');
           var googleLayer = new OpenLayers.Layer.Google("Google", { 'sphericalMercator': true,
           'maxExtent': new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34)});
           var wmsLayer = new OpenLayers.Layer.WMS("OpenLayers WMS",                             "http://vmap0.tiles.osgeo.org/wms/vmap0", { layers: 'basic' });
           var yahooLayer = new OpenLayers.Layer.Yahoo("Yahoo");
           map.addLayers([googleLayer, wmsLayer, yahooLayer]);
           map.setCenter(new OpenLayers.LonLat(75.5, 8.5), 5);
           map.addControl(new OpenLayers.Control.LayerSwitcher());
           map.addControl(new OpenLayers.Control.MousePosition());
        }
    </script>


texts: QGIS: Edits in attribute table are not updated in Postgis database. Why?

texts:
As the title say. I edit values in the attribute tables of my QGIS project but there is no update to the database. I refreshed in QGIS as well as the database.

What else do I need to do?


texts: Using wildcards for column headers in QGIS Field Calculator?

texts:
I know you can use wildcards for fields like so (in the Field Calculator):

case when "column" ILIKE '%example%' then 1
else 0
end

Can you do something similar for the column headers themselves?

I ask this because I'm joining several layers (using the Joins properties) and it would make it easier for me to edit one set of filters instead of editing several filters since the column names will have to include the joined layer name.

I'm using QGIS 2.2.


texts: How to take advantage of the values of the image properties of Image Collections in GEE?

texts:
So, this is my code.:
var dataset = ee.ImageCollection('USDA/NASS/CDL')
                  .filter(ee.Filter.date('2018-01-01', '2019-12-31'))
                  .first();

var cropLandcover = dataset.select('cropland');
var cultivated = dataset.select('cultivated');
Map.setCenter(-100.55, 40.71, 4);

Map.addLayer(cropLandcover, {}, 'Crop Landcover');
Map.addLayer(cultivated, {}, 'Cultivation');

I want to use one of the Image Properties, cropland_class_values, so it only shows the parts of the image that have a cropland_class_values value of 1, so it only shows crop cover for corn, which has a cropland_class_values value of 1.
Here is the Cropland Image Collection from the catalog: https://developers.google.com/earth-engine/datasets/catalog/USDA_NASS_CDL


texts: Leaflet if statement for mouseout function

texts:
I have a problem. I would like to set a separate mouseout conditions for one of the geoJson.property value.

The sitoation looks as follows:

I set function for

function style2(feature) {
if (feature.properties.Post_Survey_Home_Count === "")
    return{
    fillColor: '#ffff99',
    fillOpacity: '#ffff99',
    weight: 0.5
    }
    else
return {
    fillColor: getColor2(feature.properties.Type),
    fillOpacity: 1,
    weight: 2,
    opacity: 0.75,
    //color: '#ffffff',
    //dashArray: '3'
};
}

and next in terms of these statements I would like to set a separate mouseout function for this one geoJSON property:

So that I tried:

   layer.on({
                mouseover: function (e) {
                    this.setStyle({
                        'fillColor': '#663300',
                        'weight':11
                    });
                },
                mouseout: function (e) {
                    if (feature.properties.Post_Survey_Home_Count === ""){
                    this.setStyle({
                        'fillColor':'#ffff99',
                        'weight': 0.5
                    }) else {
                    this.setStyle({
                        'fillColor': getColor2(feature.properties.Type),
                        'weight': 2
                    });
                }
                }

But my map is gone.
Console says about unexpected tokens, when I fixed them the console says: Unepxpected "else". Does it mean, that if statements are not available for the mouseover and mouseout function?

1 - before mouseover
2 - after mouseover - the view back to typical for a whole layer instead of this one property


texts:
I have a raster map without hillshades and want to add Hillshades derived from the SRTM90 dataset.
I am looking for a tool, that can handle that as automatic as possible?
I don't want to bother with too much with Data-Download, -Conversion and -Processing. I am looking for a tool, that can be handled by an end-user!


texts:
I'm trying to get some GRASS functionality in a Python Script, specification r.series, but I'm struggling to import GRASS. I'm trying to uses the method outline on the GRASS Wiki, Working with GRASS without starting it explicitly, by using this code script.setup module but with these modification:

import os
import sys
import subprocess

# define GRASS Database
gisdb = "D:\GRASS_DB"

# specify (existing) Location and Mapset
location = "LearnGRASS"
mapset = "PERMANENT"

# path to the GRASS GIS launch script
grass7bin = r'C:\Program Files\GRASS GIS 7.4.4\grass74.bat'

# query GRASS GIS itself for its GISBASE
startcmd = [grass7bin, '--config', 'path']
try:
    p = subprocess.Popen(startcmd, shell=False,
                     stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
except OSError as error:
    sys.exit("ERROR: Cannot find GRASS GIS start script"
         " {cmd}: {error}".format(cmd=startcmd[0], error=error))
if p.returncode != 0:
    sys.exit("ERROR: Issues running GRASS GIS start script"
         " {cmd}: {error}".format(cmd=' '.join(startcmd), error=err))
gisbase = out.strip(os.linesep)

# set GISBASE environment variable
os.environ['GISBASE'] = gisbase

# define GRASS-Python environment
grass_pydir = os.path.join(gisbase, "etc", "python")
sys.path.append(grass_pydir)

# import (some) GRASS Python bindings
import grass.script as gscript
import grass.script.setup as gsetup

# launch session
rcfile = gsetup.init(gisbase, gisdb, location, mapset)

# example calls
gscript.message('Current GRASS GIS 7 environment:')
print(gscript.gisenv())

gscript.message('Available raster maps:')
for rast in gscript.list_strings(type='raster'):
    print(rast)

gscript.message('Available vector maps:')
for vect in gscript.list_strings(type='vector'):
    print(vect)

# delete the rcfile
os.remove(rcfile)

However, I am getting this error message:

  Traceback (most recent call last):
    File "D:/Training/LearnGRASS/LearnGRASS_7.py", line 28, in
      gisbase = out.strip(os.linesep)
  TypeError: a bytes-like object is required, not 'str'

It apppead that the offend line was this one:

gisbase = out.strip(os.linesep)

where the script is trying to remove the endline characters from out, so I modified the line to:

gisbaseStr = str(out).replace("\\n'", "").replace("b'", "")
gisbase = str.encode(gisbaseStr)

However, that return the error message:

  Traceback (most recent call last):
    File "D:/Training/LearnGRASS/LearnGRASS_7.py", line 34, in
      os.environ['GISBASE'] = gisbase
    File "C:\Users\bradl\Anaconda2\envs\GeoSpatial_3.6\lib\os.py", line 674, in setitem
      value = self.encodevalue(value)
    File "C:\Users\bradl\Anaconda2\envs\GeoSpatial_3.6\lib\os.py", line 730, in check_str
      raise TypeError("str expected, not %s" % type(value).name)
  TypeError: str expected, not bytes

But leaving it as a string kick up different Error Message:

  Traceback (most recent call last):   File "C:\Program Files\GRASS GIS 7.4.4\etc\python\grass\script\core.py", line 39, in
      import builtin ModuleNotFoundError: No module named 'builtin'

  During handling of the above exception, another exception occurred:

  Traceback (most recent call last):   File
  "D:/Training/LearnGRASS/LearnGRASS_7.py", line 41, in
      import grass.script as gscript   File "C:\Program Files\GRASS GIS 7.4.4\etc\python\grass\script__init__.py", line 5, in
      from .core   import *   File "C:\Program Files\GRASS GIS 7.4.4\etc\python\grass\script\core.py", line 44, in
      from os import environb as environ ImportError: cannot import name 'environb'


texts:
I have a MultiLinesString layer in a geopackage that I want to import to a Garmin GPS as a track. I exported the layer from QGIS as a .gdb. However, it's too long to import to the GPS (~2500 points, limit is 500 points.) I try to simplify it using gpsbabel:

  gpsbabel -i gdb -f ContDiv.gdb -x simplify,count=475 -o gpx -F ContDiv3.gpx

The output contains all the original points, all 2500 of them. The first sentence of the simplify docs say that the simplify option simplifies routes and tracks, but then the rest of the docs only says routes. Is my problem that gpsbabel only simplifies routes, not tracks?

Is there an alternate pathway from geopkg to gpx or gdp that will simplify the line string into a Garmin track, perhaps using QGIS? I know the processing toolbox has a simplify operation, but that is based on distance error, so it's a matter of hunt and peck to find a parameter that reduces the line string to 500 points.


texts: Drawing Polyline with user input points in QGIS using Python?

texts:
I have an empty layer in QGIS in that i just want to draw a polyline by selecting points on layer using Python.

For Example if i select a start point and end point on a layer , i want to draw a line connecting to those points using python code.
i have gone through all related references but i didn't get it.

what i have tried is :

from qgis.gui import QgsMapToolEmitPoint
def display_point(point, mouse_button):
  coords = "Map Coordinates: {:.4f}, {:.4f}".format(point.x(), point.y())
  print coords
  layer = iface.activeLayer()
  feats = [ feat for feat in layer.getFeatures() ]
  geo_pt = QgsGeometry.fromPoint(QgsPoint(point.x(), point.y()))
  id = -1
  for feat in feats:
    if geo_pt.within(feat.geometry()):
      id = feat.id()
      break
  if id != -1:
    print feats[id].attribute('name')
  else:
    print "no feature selected"
canvas = iface.mapCanvas()
pointTool = QgsMapToolEmitPoint(canvas)
pointTool.canvasClicked.connect(display_point)
canvas.setMapTool(pointTool)


texts: What is causing sr.lock on tables after looping?

texts:
The following code puts an sr.lock on each table in my List.   What’s causing my tables to lock up?

I understand my script is reading or opening my tables temporarily to perform the “MakeRouteEventLayer_lr” Geoprocessing but since I’m not necessarily inserting a cursor I should not have to delete references to a table right?

As a result of my tables locking up while running the “MakeRouteEventLayer_lr()”,  I am not able to run the second Geoprocessing tool “arcpy.AddXY_management()”.

Currently, in order to accomplish my task I have to run each of the Geoprocessing as a standalone process and even then the tables are locked once they finish processing.

I would greatly appreciate it if someone can help me isolate this problem. By the way my source data is in a file geodatabase. I've ran the code below in pythonWin and ArcMap with the same results.

# First Geoprocessing

Table_List = arcpy.ListTables()
for t_list in Table_List:
    arcpy.MakeRouteEventLayer_lr(PODS_ROUTEM, "ROUTE_CODE", t_list, "Route_Code POINT measure", t_list+'_EVENT', "", "NO_ERROR_FIELD", "NO_ANGLE_FIELD", "NORMAL", "ANGLE", "LEFT", "POINT")
    arcpy.CopyFeatures_management(t_list+'_EVENT',t_list+'_EVENT')

# Second geoprocessing

fc_list = arcpy.ListFeatureClasses()
for fc in fc_list:
    arcpy.AddXY_management(fc)
    arcpy.CopyFeatures_management(fc,fc)


texts:
I am trying to append one polygon shapefile to an empty polygon feature class. The append works without any errors that I can use to troubleshoot. Unfortunately, out of 59 polygons in the shapefile, only 52 are added to the empty feature class.

I'm not sure what steps to take to troubleshoot this issue. Does anyone have any suggestions on what could be causing the problem?

Edit: Steps taken so far;
I tried appending just the 7 missing polygons. The Append tool runs fine and reports no errors..but does not add the 7 missing polygons. I tried exporting those polygons to a new shapefile and appending that shapefile, also to no luck.

I ran Check Geometry and found no geometry errors. All 7 parcels have only the default ArcMap fields (objectID, shape, Shape_length, and Shape_area) plus one ID field which is just 0 for every polygon including the ones that appended successfully.


texts: Display freight volume by city on a geographic map

texts:
I have been requested to show on a geographic map the volume (in meters cubed m3) of freight (parts) that ship to various cities (Canada/USA).
The idea is to show which cities receive the most freight by volume.

Presently the only mapping tool I have accessible is Microsoft Mappoint 2002.

Is there another webbased tool that I could use that does a good job of showing volumes relative to each other?
Mappoint has a feature to show sized pie charts  - is the best type of method of showing total volume by city along with the individual contributors?

Update1:

I am not trying to show the flow of parts - only the size of origin and destination.  For example does Toronto, Ontario have 100m3 of parts and Buffalo, New York only 30m3 of parts?


texts:
I was able to use the Points to Path plugin to create a simple line with a label in the middle:

Here is the very simple text file I used:
ID;GEOM
"BL2990";"LINESTRING(490206.73 214102.673, 490236.94 214134.576)"
The problem is I can't recreate this again? I get the following error:

Dependent field could not be populated.

How do I fix this error?
This might provide some clues but I can't figure it out:
https://issues.qgis.org/issues/21055
Perhaps I have to create a layer first.
I like this approach because it is one file and I don't have a lot of steps like importing a text delimited file and then have to connect them and such.


texts: How to Classify With Multiple Ancillary Datasets using QGIS?

texts:
I've tried using the classifiers in ArcMap which only supports one ancillary dataset and I tried using Orfeo in QGIS 2.18, though the Orfeo Image Classification tool does not accept any ancillary raster datasets. I can create a training model using multiple bands and images in Orfeo but Image Classification only accepts one raster file so I can't feed it the extra ancillary data I used during the training.

Besides the original satellite image, I want to add in NDVI, NDWI, and DEM's to help the classifier just so you guys have an idea of what I'm trying to do. I haven't come across something like this besides creating a completely custom classifier.

How do I Classify With Multiple Ancillary Datasets using QGIS?

EDIT: I've figured it out. I just make a new composite raster image which includes the ancillary data as additional bands or create a virtual raster catalog. This way I can use the new composite raster with classifiers in any Remote Sensing software which often accepts only 1 raster file at a time. I've already tested this with my low spectral resolution data and it helped my classification output by adding NDVI and a NDWI as extra bands.


texts: How to get a date string to run in Time Manager?

texts:
I have in my Attribute Table two columns, one with dates and one with times (running in increments of one hour), so each feature has a basically a date and hour.

The example above shows multiple features at one point in time, but there are thousands of features and about 30 hours/timesteps.

I'm trying to get QGIS' Time Manager to read both into itself.

I tried creating a new column of data. I created a 22-character wide string and populate the column with

concat("DATE",' ', "TIME")

and it outputs the following

The second column I used

concat("DATE",'T', "TIME")

just to try this way as well.

So the output looks okay and I try to load the layer into QGIS, choosing one of my new columns. It defaults to "SECONDS FROM EPOCH (float)" and the slider starts at the beginning of UNIX time, 1970-01-01.

When I try to force Time Manager to '%Y-%m-%d %H:%M:%S' or '%Y-%m-%dT%H:%M:%S', it fails and gives me the following error.

I saw this question from several months ago that seems to be similar to my problem but it doesn't look like they got a solution. Wondering if anyone can help me? I thought strings would work but perhaps I'm missing something.

EDIT: Here are the fields I have and their properties. Is there a difference between 'String' and 'string'? I created all the original fields in Fortran and I notice the new fields I made in QGIS are labelled as 'string'.

Interestingly, if I set Time Manager to treat the data as "SECONDS FROM EPOCH (float)" (i.e. manually setting it, rather than letting it figure it out), I get the following error

Some data to reproduce the issue can be found here, I'll leave it for about a day.


texts: How to generate TFW files for a list of images?

texts:
I have a few digital images that are not georeferenced. E.g.:

Img01.tif
Img02.tif
...
Img99.tif

On the other hand, I have a text file, which contains georeferencing information for these images. E.g., I have a file c:\tfws\list.txt which contains the following information:

Img01.tfw 0.25 0 0 - 0.25 456258,125 4569852,125
Img02.tfw 0.25 0 0 - 0.25 456586,125 4570001,125
Img03.tfw 0.25 0 0 - 0.25 456952,125 4570300,125
...
Img99.tfw 0.25 0 0 - 0.25 458412,125 4575123,125

I need to generate a TFW file for each image in the text file. E.g. create a file named C:\tfws\img??.tfw containing the following information:

0.25
0
0
-0.25
456258,125
4569852,125

Does anyone know of a tool able to generate these TFW files automatically?


texts: Creating multiple routes at same time with different stop locations using ArcGIS Network Analyst?

texts:
I have a feature class that contains about 500 cross-streets locations and another feature class that contains about 1800 locations.  I need to create a route from each cross-street to 1-6 locations.  Is there a way that I can set up Network Analyst or ModelBuilder to solve for multiple routes at the same time?  For example, cross-street A route to location 1, 2, 3; cross-street B route to location 4, 5, 6 etc.

There is an established relationship between cross streets and locations, and there is an attribute that defines that relationship. However the relationship is one to many: one cross-street to multiple locations that need to be defined by separate routes to each rather than just one route with multiple stops.

I was able to model the network but am having a hard time figuring out the iterate.  I keep generating networks for all points, not just the ones I specifically need: point 1 with point A; point 2 with point B, etc.  Here's a screen shot of the model builder showing routes I was able to create between Bell2 and OID 533; and Bell2 and OID 534.  Also on a side note, how can I change the attribute fields "FirstStopID" and "LastStopID" to recognize names of the attributes of the stops--e.g. Bell2 and 533.


texts: Which python version for standalone script with Qgis module

texts:
I'm writing a standalone python script, and I need to use Qgis module.
For information i'm working on a Mac and I already have Qgis 2.18 installed.
So in list of module importation, I wrote this :

import qgis

But when I run my script, it returns me :

No module named 'qgis'

I search where it comes from, and apparently, it's because I have many python versions installed. And I dont use the one who is connected to Qgis.

I check my python install with this command which -a python python3and it returns :

/opt/local/bin/python
/usr/bin/python
/Library/Frameworks/Python.framework/Versions/3.6/bin/python3
/usr/local/bin/python3

According to python --version I'm using Python 3.4.7

Somebody can help me to use the good one ?


texts: convert from PCS: Palestine 1923 Palestine Grid to WGS84 coordinate system

texts: I have the coordinate system of points in GIS in PCS: Palestine 1923 Palestine Grid, is it possible to convert the coordinate system to WGS84 coordinate system?


texts: Tool for adding buildings to OpenStreetMaps by tracing their countours on satellite maps

texts:
There is an OSM map of some city, in which several buildings are missing (see the OSM map and the Yandex map of the same place below - there are lots of spaces, which look empty in the OSM map, but have buildings in the Yandex map).

I want to add missing buildings to the OSM map using a tool, which

displays the OSM objects over a satellite map (such that I can see what objects are present on the satellite map, but not in the OSM map),
allows to add new objects (ways) to the OSM map by clicking points on the satellite map and connecting them to a polygon and
allows to save my changes (coordinates of newly added buildings) locally in a file.

What tool are there, which support these functions?


texts:
I am currently working with QGIS 2.0 in an Ubuntu environment. I have a quite large project with many different layers and layer-groups. By default the layernames are in bold fonts and the groups are not. I guess it would be more convenient if it would be the other way around.

Is there any possibility to change it?

The setting in "general" -> "bold group box titles" is not the thing I want, as it makes the font within the settings menu for groups bold.


texts: Computing the line of equal distance among several polygons in R

texts:
I have a simple GIS problem I was not able to find the solution on the internet.
My aim is quite simple: I have a polygons layer that contains several entities, and I would like to compute a new layer that would contain a line that separates the space surrounding each entity based on an equal distance from the closest entity.
To illustrate the question, I have drawn a schematic that sums up the problem.

Figure A represents a given shapefile that contains four polygons.
Figure B includes the second layers (dashed blue line) that separate each entity based on equal distance from each other.
Ideally, I am looking for an answer on R or QGIS, as they are the two software I usually use.


texts: Using well survey data to pinpoint wells from BLM section division

texts: I feel like there has to be a good way to do this in ArcGIS but can't seem to find it.  I have well survey's that pinpoint well locations using the FSL (From southern line of section)/FWL (From western line of section) convention.  These descriptions are essentially coordinate pairs but I can't seem to find a good way to accurately use them in ArcGIS.  I can use the measure tool and get a good estimate, but I would love to have a way to specify the BLM section and then enter the FSL/FWL and pinpoint the well. Does anyone have any ideas of how to do this?


texts: How to preserve arc linetype when exporting from .dwg to ArcGIS for Desktop?

texts:
I have a closed polyline in a .dwg that I need to export to a feature class as a polygon.

The polyline is composed of straight segments as well as arcs.

The arcs, however, get converted to many short, straight line segments when exported to shapefile.

It is required that the arc line type be preserved.

How can this be done?

I am using/have access to AutoCADMap3d, ArcGIS 10.2.2, and FME Desktop 2014 tools.


texts: How do I subdivide and edit polygons preserving topology?

texts:
I'm using QGIS. I have a layer made up of many adjacent polygons - let's call them "countries". I currently use the Atlas plugin to export each "country" polygon to a .png file.

I want to further subdivide each polygon into 3 or more sections (let's call them "states"), and to be able to export these in the same way. However, I want my original polygons ("countries") to remain intact (and hence exportable). Furthermore, if I alter my "country" boundaries, I want the "state" boundaries that are affected to also be automatically updated.

Can anyone advise me on a way I could achieve this?

Thanks.

EDIT: Here's a rough mock-up of my state/country relationships:


texts:
I want to join some fields from a .geojson file wihtout geometry to another .shp layer.

The join function is not working because I want to keep the fields in the .shp file even then the .geojson has been removed from my QGIS project. And I don't want to export and save the .shp file as well.

Therefore I'm searching for a way to permanently join the fields.


texts: Is it valid to resample a raster from 1m to 3 feet?

texts:
I know that in general it is not a good idea to up-sample. I am curious if it would still be be a "bad idea" with a minor change like this from 1m resolution to 3 feet? Going to 3 feet would not actually be up-sampling anyway. Three feet should be slightly less than 1 meter. Would you recommend against this or is this perfectly valid? I am contemplating doing this with a 1M DEM.

I found many questions about resampling but most involve obvious down sampling. This would be more like changing projection from a Raster in units of meters to a raster in units of feet and rounding to 3 feet (that is how I thought to proceed).

If there is error involved, is there any way (preferably easy) of quantifying how much has been introduced by this method?


texts: Lost geometry from QGIS mysql vector layer

texts:
I have a vector layer displaying geometry from a mysql view.

This works fine for several days/weeks until suddenly it fails to refresh and my layer reverts to a non-geometry table, losing all symbols, labels etc.

Closing the project without saving and reverting does not revert back to the previous working state. I am able to add a new vector layer to the same view that displays correctly but without my custom symbols.

Is there any way to fix my layer (re-establish the geometry field) without having to start from scratch?

Edit: QGIS 2.18.16, MySQL 5.6.14


texts: Is there a GeoDjango tutorial for GeoRSSFeeds?

texts:
I'm going through GeoDjango basics right now. Ideally, I'd like to publish a GeoRSS from a set of data I have that is already in Django's database. I know how to get a regular Feed working, but when I try GeoRSSFeed, I'm running into problems.  One of my main problems is:

Exception Value: __init__() takes at least 4 arguments (1 given)

In reference to my urls.py file when I try to make my Feed to GeoRSSFeed. Is there no tutorial for getting GeoRSSFeeds to work? I've seen this but it really isn't helping me.


texts: Openlayers 3 - click under geolocation layer

texts:
I am using geolocation feature of Openlayers 3.13, as in the examples. (http://openlayers.org/en/v3.13.0/examples/geolocation.html)

However I also have queryable point layers, and when I track the position, I cannot click under the "accuracy feature".

Do you know any way of avoiding this problem?

The best would be to not detect any feature at all for the accuracy (even not change the cursor).

I tried z index and render order, but they did not work.

I ended up removing accuracy, but I would like to find a solution.


texts: Does ArcGIS for Server support Windows OS In Place upgrades?

texts: I have ArcGIS for Server running on Windows Server 2008. The IT department wants to upgrade the server to Windows Server 2012. Does ESRI officially support this migration path? I know it was possible before but with the prevalence of VMs nowadays I'm not sure it's still supported. I can't find anything official on their website.


texts: Warning message in raster::rotate(): "this does not look like an appropriate object for this function"

texts:
I'm attempting to rotate a stacked raster using the raster package in R. The raw stacked raster looks like this:

When I attempt to rotate the raster to change the longitudinal scale from 0-360 to -180-180, I get this warning message:

Warning message:
In .local(x, ...) :
  this does not look like an appropriate object for this function

However, the transformation does go forward, giving me:

What's the warning about?


texts: API or lib to convert boundary boxes or single long lat to quadtile?

texts:
I would like to store geo coordinates as OpenStreetMap quadtiles. Is there more doc somewhere than in the link above? Are there standard solutions to convert long / lat points and boundary boxes into their quadtile coordinates?
(I code in Java).


texts: Remove polygons having minimum length and minimum width less than a certain number

texts:
I am using QGIS 2.8.6 and I have a shapefile with polygons. I want to remove all the polygons which have both minimum length and minimum width less than 500m. I need to do it within a Python script. Does anyone have an idea about how to accomplish this?  Can you also explain to me how to do it with the available functionalities of QGIS?

When I tried to use Polygons from layer extent, I could find the height and width of the bounding boxes greater than or equal to 500m however these bounding boxes represent maximum length and width so I can still have polygons inside which do not satisfy my requirement.


texts: Geoserver rendering hindi labels incorrectly

texts: I have published a layer in geoserver with database. Unfortunately geoserver is not rendering the label correctly as its in database. I have included Hindi font both in sld and geoserver but I was unsuccessful. Please provide some way out to over come the above problem. Can you tell how a hindi font can be used in geoserver.


texts: QGIS : Creating point feature from Lat, Lon SQL Server 2008R2 data import

texts:
I am importing a table having no spatial columns from SQL Server 2008R2 database into QGIS 2.10.1. The table has latitude & longitude values. How can I create point features using the latitude, longitude data in the table in QGIS?
CSV import is one way but can I do it without CSV import directly from imported data?


texts:
I am using Python 3.7 and QGIS 3.4 to run a number of processes, and while I draft and debug my script I want to be able to occasionally stop, check if I want to continue, and then do so.

Usually, I would use input as so:

check = int(input("press 1 to continue"))
if check == 1:
     print("continuing")

This results in:

RuntimeError: input(): lost sys.stdin

However, this isn't working in the python console in QGIS. It seems like this might be a known bug in QGIS, but any answers I've found are pretty old and haven't worked, or are geared towards people creating plugins and actual UI.

I'm not looking for anything fancy - like I said, just using this for rough scripting and debugging on my own end.


texts: Extracting bigger raster with smaller raster in way that bigger raster has NoData where smaller raster cells are?

texts: I want to extract the raster cells from a big raster with the raster cells from a smaller raster which is overlapping with the bigger raster. Unfortunately, I cannot find a tool for that in ArcMap 10.3.1, but is there a way to do it? I do not want to subtract the values, but there should be noData in the new bigger raster where the bigger raster is overlapping with the smaller raster, like a hole in the bigger raster with the extent of the smaller one.


texts: How to add images in CartoDB infowindows

texts:
Can't seem to get the default site to display images.

I've added HTML anchor hrefs to a image column, for example:

I've dragged the image column to the top of the infowindow sidebar app and it still just displays "Non-valid picture URL". I also can't seem to figure out the last sentence of the explanation.

The explanation is:

Adding images to your infowindows
This template allows you to add an image to your infowindow.  The image URLs will need to be stored as a column in your table. Next, you will need to drag the URL field to the top of your infowindow and then select 'Header with image'.


texts: Building table about LULC historical (rows: categories; column: data by year)

texts:
I'm working with datas of Land Use and Land Cover (LULC) (annual shapefiles) in a period of 1985 to 2018. I want to build a historical LULC table (rows = category; column = year LULC) using these shapefiles. I need to build this table for many watersheds. Nowadays, I'm building this table manually.
Are there any tools to build a table by shapefile attributes in QGIS?


texts: ArcGIS Online and Sign In options greyed out in ArcMap?

texts:
I have just installed ArcGIS for Desktop 10.2.2 on my Windows 7 machine. When I got ArcMap opened up I found the "ArcGIS Online ..." and "Sign in..." options were greyed out. Meanwhile in ArcCatalog there was no GIS Server -> ArcGIS Online option.

Just wondering if there is any options/settings I should tick on to get them working?


texts: Land Use Public Notification: Invalid geometry

texts:
I am in the process of incorporating our data into the Public Notification app from ESRI, but when I load in our parcels, I get an error saying "Unable to perform operation: invalid geometry". I have checked the spatial reference, and it is the same as the sample data, and I have tried repairing the geometry, and the problem is still around.

Would I have to transfer my data into the Local Government data model?

Also, I have checked the Chrome DevTools, and all it tells me is "Error: Failed to execute query".


texts:
I have the following data of a GEOTag photo

GPSMapDatum     WGS-84

GPSLatitudeRef  N

0)              18/1

1)              57/1

2)              2147483647/89750086

GPSLongitudeRef E

0)              74/1

1)              10/1

2)              2147483647/54998752

GPSAltitudeRef

GPSAltitude         2147483647/3510369

GPSImgDirectionRef  T

GPSImgDirection     2147483647/81244832

My question is how can i convert this data to proper lat/long co-ordinates
(eg: 18.1234, -74.1234).


texts: QGIS Geometry Generator - Connecting Points with Lines - with a 'one to many' relationship

texts:
I'm working with  QGIS 3.12 and have two large point shapefiles (buildings/barns and wetlands) that I want to connect with lines via the geometry generator. I've found a function called get_feature() that will connect the points in the correct way but the get_feature() function only selects the first in a one to many relationship.
See example screenshot. The blue wetland dot is only connecting to one of the four building dots.
The script I have currently is:
make_line(    geometry(get_feature('barns', 'wID', '3809' ))  ,  geometry(get_feature('wetland_matrix', 'InputID', '3809' ))      )

I need to connect all four, not just the one!


texts: Merging overlapping thematic rasters using ArcMap

texts:
I want to merge two overlapping rasters. Raster A has 8 classes and raster B has only 1 class. I would like to produce Raster C with 9 classes, where Raster B would overwrite corresponding pixels in Raster A. I tried using the "raster merge" function in QGIS, but it does not produce what I want.
How do I do this in ArcMap?


texts: Find higher elevation (10m) perpendicular point from stream network at regular interval

texts:
I have a stream network and for that network, I want to find the points at regular intervals which lie at 10m distance perpendicular from the stream.The final goal is to find the buffer of that 10m elevation around the stream.

I have tried hub distance, elevation, point sampling plugins, etc.
One of my approaches was to convert the DEM to points then convert the stream network to points. Then I used point sampling to attach the stream points their respective elevation information. Now I want to make a line between the stream points and DEM points based on higher elevation and perpendicular distance.

Is there any plugin or code which can help to achieve? Also, I think that there must be a better way to find higher elevation info from Vector (Stream) and Raster (DEM) layer and the buffer if you can point to any such plugin/code/algorithm that will be a great help. I need to automate the process as well.

Sample image for reference:-


texts:
I'm doing a project using the combination of GeoDjango, PostGIS, and OpenLayers

So far, I've imported OSM data with osm2pqsql, done settings and I knew the data is constructed by planet_osm_line, planet_osm_polygon, etc.

Now, I want to query all the villages boundaries under certain city and render with OpenLayers, I know I have to get the GeoJSON from the database and pass it to the template, but I have no idea how to make the queries using GeoManager, can somebody give me some instructions

Note: English is not my mother language, I try my best to explain my question


texts: Do somebody know how factors as "age of data", "scale" and others influence the precision of GIS data?

texts:
Several factors may cause imprecision in datasets. As example of such sources, we should list: "age of data", "formatting", "scale", the technique employed to create the data and others. Each of them disturbs the precision in a different way and magnitude.

For instance, the scale in which the data are produced cause a well known imprecision pattern and it also has standards like that proposed by the United States Geological Survey:

"requirements for meeting horizontal accuracy as 90 per cent of all measurable points must be within 1/30th of an inch for maps at a scale of 1:20,000 or larger, and 1/50th of an inch for maps at scales smaller than 1:20,000."

Accuracy Standards for Various Scale Maps

1:1,200 ± 3.33 feet
1:2,400 ± 6.67 feet

1:4,800 ± 13.33 feet

1:10,000 ± 27.78 feet

1:12,000 ± 33.33 feet

1:24,000 ± 40.00 feet

1:63,360 ± 105.60 feet

1:100,000 ± 166.67 feet

However, in my researches I have found no other indicators of the magnitude of errors associated with the other factors and sometimes it seems to me that nobody has done this yet.
My goal is to build a probabilistic model for GIS positional errors and use it to create probabilistic versions of spatial operations like spatial joins, for example.

Do somebody know some works that have quantified GIS erros such that it can help me achieve that goal?


texts: Debugging TypeError: SearchCursor() got unexpected keyword argument in ArcPy script?

texts:
I tried to use SearchCursor(points, fields=fieldname) but I get this error:

  TypeError: SearchCursor() got an unexpected keyword argument 'fields'

I'm not sure why this doesn't work.

    # Local variables:
    points =  arcpy.MakeFeatureLayer_management("C:\\alltogether.shp")
    homebase_shp =  arcpy.MakeFeatureLayer_management("C:\\homebase.shp")

    fieldname = "Device"
    Device = arcpy.AddFieldDelimiters(points, fieldname)
    delimitedfield = arcpy.AddFieldDelimiters(points, fieldname)
    devices = set([row.getValue(fieldname) for row in arcpy.SearchCursor(points)])

    for i in devices:
            attr_sel = arcpy.SelectLayerByAttribute_management(points,"NEW_SELECTION", delimitedfield + " = " + str(i))
            totalpts = arcpy.GetCount_management(attr_sel)

            # Process: Select Layer By Location
            inside_sel = arcpy.SelectLayerByLocation_management(attr_sel, "WITHIN", homebase_shp, "", "SUBSET_SELECTION")
            insidepts = arcpy.GetCount_management(inside_sel)

            print "device %d has %d total points and %d inside the base" % (int(i),int(totalpts[0]),int(insidepts[0]))

a comment I saw earlier suggested doing it with fields=

  Since you are only interested in one field, for better performance I would specify the >optional fields argument, e.g. myList = set([row.getValue(fldName) for row in >arcpy.SearchCursor(fcName, fields=fldName)]) – blah238 Jan 29 '12 at 1:03

Select distinct values from a single column of an attribute table (or layer)


texts: Making all layers in an SVG visible when exporting from QGIS?

texts:
My QGIS project has multiple layers. I use a U.S. county map as my base, and higher-level layers have lines. When I use the simpleSVG plugin to export my view to an SVG, then open the file in a web browser, I only see the basemap. When I use the browser's web inspector to check the SVG file in detail, I see the SVG elements for the lines, but they appear invisible, not colored in.

What do I need to do to make QGIS export all layers with visibility?

Link to the SVG code I get when exporting only the lines: https://pastebin.com/36pnwRmZ. They appear invisible when I open the SVG file in browsers.


texts: ArcSde Oracle Maximum number of grids per feature (8000)

texts:
I have a problem with the import of data (lines) to ArcSde Oracle.

2011-07-14 12:10:04| 22.9| 0.0|ERROR |Error while executing stream insert operation for table 'CHEMN'.:SDE_ERROR_CODE(-51) Underlying DBMS error. SDE_EXT_ERROR(29875) SDE_ERROR_MSG1(ORA-29875: échec d'exécution de la routine ODCIINDEXINSERT
ORA-20092: Maximum number of grids per feature (8000) exceeded.
ORA-06512: à "SDE.ST_DOMAIN_METHODS", ligne 1982

Have you a solution ?

Thanks


texts: Can I use in a postgis db within the same table a combination of CRS?

texts:
I'm experimenting around with my postgis2.0 database. I want to create a basic table to store all my point data and in the process,I realized that the coordinates are in different CRSs. So naturally I'm wondering now if the  the_geom column can contain a variety of CRSs - or Im confined to a single CRS per table.

To visualize it a bit, Is the following possible?

point_id | the_geom      |
1        | (geom_string1) <- CRS EPSG:4326
2        | (geom_string2) <- CRS EPSG:2100


texts:
I want to create a big image of this area. But with a bigger zoom level. How can I achieve this?

Currently I would zoom in, make a screenshot and put those screenshots together with a image editor. But this can't be the fastest solution for this task.


texts: A saved QGIS project doesn't show any layers

texts:
I was working with a lot of layers on a QGIS project. I saved it and when I opened the file it doesn't show any layers. The .qgis is not empty, weighing 1,7Kb (I have other with 500Kb working well) but when I open the project it only shows the names of groups I made to manage the layers on Panel Layer.

I have try to use older versions of the project and they dont work either.

How can i fix it?


texts: Turning off all layers in MXD using ArcPy?

texts:
I tried this standalone script to turn off all layers in the MXD:

print "This script turns off the following layers:"
    try:
        arcpy
    except NameError:
        import arcpy
    mxd = arcpy.mapping.MapDocument(r"PATH")
    for lyr in arcpy.mapping.ListLayers(mxd):
        print lyr.name
        lyr.visible = False
    mxd.save

It is adapted from the script mentioned for the Python window here:

https://geonet.esri.com/thread/16415

The standalone one I created for IDLE from this, above, did not affect the visibility status of any layers when I opened the MXD afterward to check. This is strange because running the same script within the Python window turned off all the layers for me.

Any clues?


texts: Identifying and quantifying gullies and hags from aerial imagery in QGIS

texts:
We are involved in an extensive peatland restoration project in the Shetland Islands. We provide high resolution and accurate digital surface models and orthomosaics for initial analysis and on-going monitoring. One of the most arduous tasks on a new site is identifying the gullies and hags (think of a one-sided gulley); this is done manually by drawing shapefile line strings using hillshade and slope rasters overlaid on the orthomosaic (see attached). What we would like to do is semi automate this process, identify and create shapefiles of gullies and hags more than 0.5m deep and 2m long.
The example image shows a small section of hillshade extracted from the 0.5cm cell size DSM, and a slope analysis using the raster calculator to display slopes over 25º. The hillshade is set to 50% transparency, and the slope is set to display white as no data and overlaid on the 0.25cm orthomosaic. It's very easy to see the gullies this way, but with the average site being between 50 and 100 hectares it's a laborious process!

The r.geomorphon tool is getting towards the results I want. The image below shows the DSM resampled to 0.1m cell size and processed using the following parameters in geomorphon: outer search radius: 35, inner search radius: 9, flatness threshold: 3.7, flatness distance: 15.

Further processing with the raster calculator isolates the gullies:

Now I need to simplify and thin the raster then convert it to a vector, however r.thin fails.
How can I convert the raster to a vector?


texts: Combining polylines if they meet specific conditions?

texts:
I want to combine the polylines in one shp if they are connecting, having the same names and same values. The same names and same values are based on two fields of the table, as shown below. Is there any existing built-in tool I can use?

I am using ArcGIS 10.3.1 Advanced License.


texts: How to find WFS layer geometry type, when the layer is empty (0 features)?

texts: Is it possible to find the geometry type of a layer in a WFS, using OGR, when the layer is empty (0 features)?


texts: Calculating loss of value from a wind turbine on surrounding properties depending on distance?

texts:
I have to evaluate where installing a turbine would be most efficient, depending on wind speed and price of properties surrounding each cell in a raster. The cost on value of properties is computed using some estimation based on the distance of the turbine to nearby properties.

I know the price of each property in the area. I have to find a way to calculate how much each turbine would cost, in terms of loss value, on the surrounding properties.

I thought of doing this by making buffers for each given distance, joining the resulting layer with the layer containing the price of properties, and then making a buffer reducing the buffer containing the value of properties to the cell (therefore, in the same, negative amount of the first buffer). However, the value of properties does not appear in the last buffer layer:

(example) Cell => 600 m buffer => join buffer with layer containing value of properties => buffer - 600 m => cell

Is there a way to keep information from the original cells on the buffers created?

Is there an alternative way to calculate the loss of value caused by a cell on the surrounding cells?


texts: I am building a plugin in python for Qgis. I want the users to be able to start typing a text and that auto compellation to pop up based on values in a certain field I will define. Is this possible and how? Is there an appropriate widget for this purpose?


texts:
I'm using QGIS 2.12 and trying to load the following WMTS layer
http://geoappext.nrcan.gc.ca/arcgis/rest/services/BaseMaps/CBME_CBCE_HS_RO_3978/MapServer/WMTS/1.0.0/WMTSCapabilities.xml

from http://geoappext.nrcan.gc.ca/arcgis/rest/services/BaseMaps/CBME_CBCE_HS_RO_3978/MapServer

When I add a new WMS layer in QGIS I can connect to the server and add the desired layer but the layer doesn't load. I'm stuck with a blank page.


texts: Topology checker does not show error for one polygon fully within another polygon as overlap

texts: Its about QGIS. In the attached screenshot, the topology checker does show error for the overlap for the two large polygons (1 and 2); but it does not show an error for the smallest (3) and the left polygon (1). I have tried both 'must not overlap' and 'must not duplicate' and still 1 and 3 don't give an error.


texts: Duplicate Layers When Using Multiple Data Frames

texts:
I regularly work with on an .mxd with multiple data frames, showing updates on one Area of Interest shapefile per data frame. This has always worked without a hitch to show updates on 2, 3, or 4 non-contiguous extents side-by-side.

Recently, I added a data frame/AOI that is contiguous to an existing one, but I still want to show them on individual data frames to provide discrete updates. Without me doing anything, I find the shapefile from one unique data frame also showing up in the adjacent data frame. When I toggle off visibility of the shapefile in the neighbor's data frame, after a click or two, it also turns off in it's native data frame. If I "remove" the alien shapefile from it's neighbor's data frame, I get a red flag error that the file is now missing in the native data frame.

Any ideas on a workaround here?


texts: Launching QGis 1.8.0 Lisboa on Windows 8 getting an error "the procedure entry point could not be located in the dynamic link library"

texts:
Possible Duplicate:
  How to fix QGIS error “Entry Point could not be located”?

I have a problem with QGIS 1.8.0 Lisboa installed on a desktop with Windows 8 Pro.
When I launch qgis.bat the message "The procedure entry point ?minimum@QAbstractSlider@@QBEHXZ could not be located in the dynamic link library C:\Progra~2\Quantum GIS Lisboa\bin\phonon4.dll" appears.

I would like to precise that when I installed QGIS 1.8 for the first time it ran correctly.

I think that some Windows update or new software installations caused the problem.

If I try to re-install QGIS nothing changes.

Is there a possibility to restore the situation?


texts: I know it's possible to make a trade area (service area) in ArcGIS. But Is it possible to make a trade area with a network data set in Grass GIS or QGIs?


texts: Calculating population density from urban raster reclasiffy

texts:
I am using ArcMap

Out of a raster, I reclassified it into only having the urban areas.

I ran this through zonal statistics and got the area of the urban area (in meters, which I presume are in square meters because it is area, correct?).
After getting the zonal stats, I joined into a census tract file which has immigrant population data (ex. dominicans, europeans,etc.).
I want to calculate the population density by square kilometer of these urban areas (1000 person/km2). What I did was taking the area field (results of the zonal stats) and dividing the presumed square meters into 1,000,000 to get the square kilometers.
After which, I divided the population group by the square km of that tract.

My question is, did I did this correctly, did I missed a step or are the results of the zonal stats not square meters?


texts: Using "snap pour point" to find a watershed of specific area starting from a point

texts:
I am trying to find everything downhill of a certain point in order to create a watershed of a specific area.

For example, say I have a point at (x,y), and I'm interested in creating a watershed that is 50 km^2 starting from this point.

One method that I have used of achieving this is placing random pour points downstream and using the watershed analysis tool to create watersheds from (x,y) to the new pour points until I find a point that gives me a watershed of 50 km^2.

Since I need to do this process ~250 times, I'm greatly interested in finding a more efficient method. One idea is using the "snap pour point" tool from the original (x,y) point to find another point that has the highest accumulation within 50 km^2. Unfortunately, the snap distance given in ArcGis is (I assumed) a radius.

Does anyone have any tips that might help me achieve what I'm trying to do?


texts: Connect point to closest feature in ArcGIS 10.2

texts:
I have some points and some features. These points fall outside the features. I want to connect each point with a line to the closest feature. How can this be done using ArcGIS 10.2?

I have checked around toolbox but didn't manage to find the right tool for it.


texts: How to configure a WMS layer to display with Google maps

texts:
I have the following code to diplay a WMS layer (in Geoserver) with an OL basemap:

var extent3857 = [-6016334, -2369948, -5843698, -2271608];

var projection3857 = new ol.proj.Projection({
code: 'EPSG:3857',
extent: extent3857,
units: 'm'
});

var baselayer = new ol.layer.Tile({
source: new ol.source.OSM()
});

var wms = new ol.layer.Image({
    extent: extent3857,
    source: new ol.source.ImageWMS({
        url: 'http://localhost:8080/geoserver/Prada/wms?',
        params: {'LAYERS': 'Prada:Pradaria', 'STYLE': 'grass'},
        serverType: 'geoserver'
        })
    });

var map = new ol.Map({
    layers: [baselayer, wms],
    target: 'map',
    view: new ol.View({
        center: [-5930030, -2316645],
        zoom: 10
    })
});

And it works fine,
the thing is I want to change/replace the basemap for Googlemap.
How can I configure the WMS or put the googlemaps as basemap?.
There's multiple non-update info sadly (openlayers-2, googlemap-v2, etc.).


texts: What are these files? Why they cannot be turned on/off?

texts:
My question is why there are some layers which cannot be turned off?

In my Google Earth project those were main maps and then other polygons/paths were, but In QGIS I am not able to turn them all off/on?

What am I doing wrong? I use QGIS 3.2 Bonn


texts: Nested iterators - batch merge from multiple outputs

texts:
I'd like to ask you some advice for nesting an iterator (from a submodel) into a master model which also include an iterator (see pic attached). The whole purpose is to 'Merge' the iteration output of the submodel with an iterated file from the master model.
The tricky part is where both iteration have to occur successively and assemble the right iterated features at the right timing before the 'Merge' operation.

Any experience in this?


texts: Trying to use Tkinter to open dialog box in ArcMap

texts:
The following script works great from the command line but crashes from the pythonaddin:

import Tkinter

import tkFileDialog
root = Tkinter.Tk()
root.withdraw()
pathtest = tkFileDialog.askopenfilename()
print pathtest

pathtest = pythonaddins.OpenDialog('Open Item', True) does not work for me because it suppresses filetypes that I need to see (pdf).

I know this is a duplicate question but I just joined the site and could not comment on other posts.


texts: Raster reprojection in QGIS leads to different results

texts:
I have a base ASTER 30m raster map from which contours have been extracted / reprojected to UTM using different processes in QGIS 2.14:

ASTER (WGS84) > save raster as using Project CRS (UTM) > Extract contours (Grey layer)
ASTER (WGS84) > Extract contours > reproject on the fly (UTM) (Yellow layer)
ASTER (WGS84) > warp raster to UTM > Extract contours (Blue layer)

I naively thought that if starting from the same geographic coordinate system (WGS84) and ending up in the same projected coordinate system (UTM) they would perfectly overlay each other but as you can see that is not the case. The different contour layers are separated from each other by a varying distance.

I'm wondering what is going on 'under the hood' to cause this... Is it a quantisation / rounding error? Does OTF reprojection use different algorithm from warp?

Is one method more accurate than the other?


texts: Why QGIS 2.12 "forgets" selection colour opacity settings?

texts:
I have few vector layers which I need to edit frequently. I'd like to be able to see through the layer when it's selected. I changed selection colour from both: Project Properties General settings and Options-->Canvas&Legend settings. For both of these I changed the opacity to 40% --> click OK --> save. It works for a while, but after I restart QGIS, the opacity is still set to 40% in both settings, but I can't see through the selection. So QGIS kind of thinks it's still showing transparent selection colour, but it actually is not. So is this a bug in QGIS or am I just missing something?

If I change the selection colour opacity from 40% to 39% and click "OK" the layer comes transparent again. I'm working with PostGIS layers, but this same problem occurs with shapefiles also.


texts: Text Background objects not overlaying properly

texts:
To complete a map of highways, I've needed to create labels with text backgrounds. I'm attempting to recreate a proper County route sign such as this:

So, I made a text background out of 3 objects: A small blue pentagon, a slightly larger yellow pentagon, and a slightly larger blue pentagon. Then I label it with the number of the route. It comes out looking like this:

You can see how the individual pentagons are not properly centered on one another. I know how to fix this manually. I can just readjust the x and y positions of the individual pieces until they line up properly. What I want to know is what is causing the individual pentagons to come out uncentered from one another. Especially since the same symbol will come out fine in one map, screwed up one way in one map and screwed up another way in another map!

Edit: I am using Maplex label engine.