texts:
I'm having a problem with heat map in QGIS. I have two shapefiles with accidents (with and without victims). Every time I try to create a heatmap with them, it gets pixilated as you can see in the image.

How can I make it smoother and not square, like in this image?

I put another picture with the options. QGIS documentation says the number of rows is connected to the number of columns. By doubling one, the other also increases. I noticed, on other maps I made, that the number of columns is greater than the number of rows. In this shapefile the opposite happens. The number of rows is greater than columns. I think that's what affects the outcome, giving this pixelated appearance. I have already tested various radius sizes (10, 30, 50, 1000, 30000 and others), but the result is always almost the same.

I'm using EPSG:4326 (WGS 84)


texts: How much information does QGIS's OpenStreetmap Plugin download?

texts:
I am using the OpenStreetmap Plugin for QGIS to download the information for my university extent.

The information as I see is downloaded in vector format divided into Points, polygons and lines...

I want to confirm whether this corresponds to all the elements of Open street map within that area with all associated tags or does it exclude something that the maps contain ?


texts: How to automatically set polygons borders to improve map readability when rendering thematic maps?

texts:
When mapping an value as a color to render a thematic map the objects borders can alter the map readability a lot. For example in this thematic map of the US counties with black borders:

There are a lot of small objects in some areas, and it makes them look darker because of the borders without respect to the thematic color. If the borders are drawn with 30% opacity it results in the following map:

This one is obviously better at communicating the data mapped on the objects. However having a static opacity rate can also degrade readability since it makes the objects boundaries less visible, which is a problem if there is only a few big objects. Moreover the change in readibility also depends on the contrast between border and fill colors.

Do automatic ways of finding the good color/opacity for layers exist? What GIS software provide method to properly handle this issue?


texts:
It seems as if Desktop 10.4 is not compatible with Kaspersky Endpoint Security 10: Whenever I want to open a saved mxd (or the application itself) the loading time is increased from the normal 10 to approximately 70 Seconds, with the splash screen. Apart from disabling Kaspersky, nothing has helped so far, e.g. adding all applications from the ArcGIS system folder to the list of save applications in Kaspersky, disabling "application control" altogether (in Kaspersky)...

Is this a known issue, soon to be adressed with an upcoming SP release? Or do I have to uninstall Kaspersky? As you can imagine, productivity suffers a great deal from this delay.


texts: Since installing QGIS 3.10.1-A Coruña I can't get Python script to work in IDLE: ImportError: No module named qgis.core

texts:
I have tried adding sys.path.append and under cmd.exe set Python Path but am still getting the dreaded ImportError: No module named qgis.core.

I set PythonPath to C:\OSGeo4W64\apps\qgis\python

I have referenced sys.path.append material on this forum and on the QGIS site and I am still getting the error.

I am using IDLE to run scripts.

Can you provide any guidance?

Here is my Script to round raster cell values:

# Description: Converts each cell value of a raster to an integer by proper rounding
# Import required modules
# Python modules

import sys
sys.path.append("C:\OSGeo4W64\apps\qgis\python")
sys.path.append("C:\OSGeo4W64\apps\Python37\lib\site-packages")
sys.path.append("C:\OSGeo4W64\apps\qgis\python\qgis")

import time
import os
# Qgis modules
from qgis.core import *
import qgis.utils
from PyQt5.QtCore import QFileInfo, QSettings

import subprocess

subprocess.check_call([
"gdal_translate", "-of", "GTiff", "-ot", "Int32",
"H:/test/demtest1.tif", "H:/test/output/int.tif"
])

#done
print("done")


texts: Set coordinate system for deed calls in COGO?

texts:
I am working with deeds that display metes and bounds, but the surveyor did not include the coordinate system used. If I am trying to enter this data into ArcMap using COGO, where in the process of entering this data should I put the coordinate system?

I am trying to add the plots of land to and existing shapefile and I want to set the data I am entering to the same coordinate system as the existing shapefile.

Any suggestions?

I am new to using COGO.


texts: Using "Difference" in Geoprocessing Tools results in Minidump

texts:
I want to use this tool on 2 vector layers but as soon as I start the process I get a minidump and I have to close QGIS.

Does the .dmp file help you in solving this? Or can I open it and look for myself? But which program do I use to open it?


texts: Does the concept of a map service exist in ArcGIS Online?

texts:
I'm new to web GIS.
I have created an ArcGIS Online account that contains:

Web Maps
Feature Layers (hosted)

It is my guess that the term feature layer is a synonym for feature service.
If that's correct, then in addition to creating feature services in ArcGIS Online, is it also possible to create map services in ArcGIS Online?

The purpose of the map service would be to query the data in the layers via HTTP.

(I suspect I've misunderstood something here. I'm finding the web GIS terminology/concepts to be difficult to parse.)


texts: ArcGIS REST API query all features inside a Polygon Geometry

texts:
I have Feature Layer A and Feature Layer B
Both layers are Poligons type. I query a feature from Layer A using que /query endpoint and get the geometry object like this
const data = {
      f: "json",
      returnGeometry: true,
      token,
      outFields: "*",
      where: `id='123456789'`,
    };

    const options = {
      method: "POST",
      headers: { "content-type": "application/x-www-form-urlencoded" },
      data: qs.stringify(data),
      url: layer_a_url + "/query",
    };

    axios(options)
    .then((res) => {
      let geometry = res.data.features[0].geometry
    })

That geometry object is something like this
{
  rings: [[]]
}

The problem I need to solve is how do I use the /query endpoint to query all the features from layer B that are inside that geometry?
I have already try something like this
const data_b = {
    f: "json",
    token,
    geometry: geometry,
    geometryType : 'esriGeometryPolygon',
    spatialRel:'esriSpatialRelWithin',
    outFields: "*",
    where: '1=1'
  };

  const options_b = {
    method: "POST",
    headers: { "content-type": "application/x-www-form-urlencoded" },
    data: qs.stringify(data_b),
    url: layer_b_url + "/query",
  };

  axios(options_b)
  .then(res => {

    res.data.features.forEach(feature=> {
      console.log(feature.attributes)
  })

But it simply returns all the features from the layer B without filtering. If I leave the where param blank or if I don't put it there, it gives me an error.


texts: I installed the OpenGeo Suite in my computer running Windows 7 Home Edition. When I installed it, I checked the box "launch OpenGeo Suite" and it ran. After installing several other software (XAMPP, Hamachi), the Suite Dashboard won't show up. As in there's no interface or taskbar icon. However, I checked the services tab in Task Manager, it was there and running. I tried changing the ports used by XAMPP, Skype, YM, and Hamachi, but the interface isn't there still. Any help is very much appreciated. Thank you.


texts: Android Tablet KML/KMZ App for polygons w/offline & GPS capabilities?

texts:
I am looking for an Android app which allows me to view and query polygons stored in KML and/or KMZ with a GPS point representing my location and offline capabilities. My KMZ files are being output from ArcGIS 10.2.2 representing some several thousand features (lakes, ponds, etc.). I would like to be able to view these features by navigating to them via GPS.

Previous poking around online has suggested that Backcountry Navigator (BCN) & Gaia GPS (Gaia) would be good places to start. I was able to load the large the large KML into BCN, but the polygons are serrated and the Object IDs are not always readable. Ideally, I could click on one of the features and look at the attributes. Gaia had issues loading such a large file, but the user interface is very agreeable.

Other posts on GIS StackExchange have suggested Locus. After downloading the free app, I found  that there were a number of issues with obtaining offline apps. For whatever reason, after trying a couple of WiFi connections, I was not able to download the tiles to my 2013 Nexus 7 (Android 4.4). Conversely, BCN & Gaia downloaded the offline maps with ease.

Does anyone have a recommendation? Bonus points if the app allows the overlay of large raster files (~200 km^2 with 1 m resolution).


texts: Misaligned .dem from Civil 3D into ArcMap

texts:
Imported 3 .dem files from a Civil 3D file. The two surfaces(topographic maps) fit correctly into my ArcMap file but the third surface, a bathymetric contour, seems to be misaligned and the scaled too big. I already tried checking the projections and units are all the same both in Civil 3D and in ArcMap. I also tried importing the files individually into a blank map but the third surface really won’t fit.

Any suggestions?

EDIT: for now I've just georeferenced the bathymetric surface into the other surfaces and it seems to be sufficient for my project but I was hoping to get a more exact answer if possible since the georeferencing is based on how precise I'm able to place the similar points and since the contours are different, it isn't really exact.


texts: What are the wider capabilities of the QGIS atlas functionality?

texts:
Given that it is already understood that the QGIS composer atlas function makes it possible to print separate pages based on a 'coverage' layer, what basic range of functions/outputs can be achieved with this tool?

The QGIS print composer atlas capability is powerful. I ask the question because even having understood the basic functionality it's really difficult to grasp the wider capabilities of this tool.

(For example, is it possible to rotate output maps, is it possible to change scales for different maps, is it possible to show and hide features, and so on?)

(I'm not convinced this is a valid question for Stackexchange - I judge that it is but I'm happy enough if others strongly disagree - the question could be re-phrased... where can I find a clear list of capabilities? I'm asking because I can't find such a list - the documentation doesn't really provide this.)


texts: What is the easiest way to let a user upload overlayers to openlayers (Geoserver/Mapserver)?

texts:
I need to let a user upload vector layers he created with QGIS or MapInfo professional to my website (using OpenLayers as client and either Geoserver or Mapserver as WMS).

What will be the easiest way for a non programmer user to draw vector layer in program (one of the two I mentioned above) and be able to upload it to my server and add it to the map?


texts: Finding latitude and longitude for any country's postal code

texts:
I have been trying to solve a problem where I have only a Postal Code and Country name, and need to know the latitude and longitude.
I have been relying on various free geolocation data dumps I've found online, but it's not common to find only postal codes and lat longs. I've also been unable to source data for many countries I am interested in (e.g., Israel).
I am exploring either paying for a large dataset of geolocation info (of which I just need postal, lat, and long, so most of the data is not useful and therefore expensive). I also will need to develop a service to ingest some custom data set.
I am also exploring using an API where I can just take a postal code and get the lat long. Paid is fine. However, all of the free APIs I've found are incomplete, and most paid APIs I've found don't expose enough sample use case or documentation for me to know if it will be a fit for my needs.
I'm pretty sure the Google API will work for me, but I simply don't have experience to know if there are better quality, or cheaper, faster, etc. solutions.
What is the most common way to take any country's postal code and get a latitude longitude for it?


texts: How to automatically refresh canvas after saving layer edits

texts: After saving a layer edit in QGIS 2.8.2 (Win7) I've to press F5 or change the map section to refresh the canvas. This is quite annoying when working with dynamic PostGIS views. Is there a way to autorefresh the canvas when saving a layer edit? Using QGIS 2.8.2 on Ubuntu everythin is working fine.


texts: How can I improve performance of gdalwarp and gdal_translate pipeline?

texts:
tl;dr: I can't get gdal_translate to use multiple cores. How to fix?

I am using gdalwarp followed by gdal_translate to process a large GeoTIFF by first cropping to a polygon cutline and outputting a virtual raster, then translating the .vrt to a .tif. I have followed suggestions from a few different answers on this site, first I split up the processes into two to enable better compression following this answer about gdalwarp, then I attempted to speed up the performance of gdal_translate following this answer about multithread support for gdal_translate. I am running this on a remote server which has GDAL v2.2.2 installed and the OS is Ubuntu 16.04.6 LTS (Xenial Xerus).

This is my code.

gdalwarp -of vrt -crop_to_cutline \
  -cutline  ${path}/counties_chesapeake_watershed.gpkg ${path}/bigraster.tif ${path}/clippedraster.vrt
gdal_translate -co compress=LZW -co NUM_THREADS=8 --config GDAL_CACHEMAX 512 \
  ${path}/clippedraster.vrt ${path}/clippedraster.tif

My issue is that I don't believe that gdal_translate is using multiple cores, though I've tried to specify this with NUM_THREADS and also to increase GDAL_CACHEMAX. This is a very large raster (~12GB, several hundred km extent at 1 m resolution) so it is running extremely slowly. Can anyone help me parallelize the compression done by gdal_translate so this will run faster?


texts: Choosing the map scale (1:nnn) for the given raster dimension

texts:
Starting points:

let define the raster image size, say: 1800x900 pixels
I want draw my own raster image map - by drawing polygons e.g. like "draw line from pixel with coordinates x100y100 to pixel x100y103"
it is easy to convert the dataset geographic-coordinates to pixel coordinates
can choose from 3 datasets, each have different scale

1:10mil
1:50mil
1:110mil

The question is:

how to choose the needed dataset (so how to choose the needed map scale)
to draw the map of the given real dimension (e.g. whole earth, or area of 2000x1000km)
with minimalizing multiple drawing of the same line e.g. to minimalize subpixel coordinates?

Probably the question is badly formulated, so trying it explain more precise:

having the earth "width" is approx. 40.000km.
e.g. for the image of the width 1800 pixels, 1 pixel is 22km.
So, drawing the world map to 1800x900 image from any data what has higher "resolution" than 22km/pixel is pointless, because will need to draw a line from a pixel coodinates x100y100 to x100y100 (e.g. from the given pixel to the same pixel, because the map data is unnecesary too much precise and after the conversion to pixels the result gives the same pixel)
so, how to determine the needed dataset scale?

Hope now it is more clear - sorry, my english is not my native...

I was read many questions here, like:

Dimensional analysis/conversion problem: Map scale to image resolution (dpi)
How to Determine to world map size in pixel given this data?
Mapserver, TMS, Scales zoom levels and how to deal with all of it
and many others - but still don't figured out an reasonable "algorithm".


texts: Converting Raster to Polygon/Shapefile using ArcMap without Spatial Analyst or 3D Analyst

texts:
I have a floating point raster file which I want as a polygon for visual purposes so I can trace/trim/etc. the polygon (unless there is a way I can trace the raster?).

Looking into other posts on here most tools recommended require spatial analyst or 3D analyst which I do not have access to.

I don't need to retain any information from the raster, just an outline of the shape itself but as it is a large complex shape it would not be easy to manually draw so is there a way to convert it to a polygon?

I am using ArcMap 10.1


texts: How to plot a full hemispheric orthographic projection in R

texts:
I want to make a plot like that using R:

All I want is the country borders as seen from outer space, no fancy colors needed. It seems I need to use the spTransform function, but it gives me a "non finite transformation detected" error. I've been told that I need to apply a cosine formula to each point in my world shapefile, and erase those that goes behind the Earth. This has 2 problems:

It looks like something the algorithm should do automatically;
I cannot just erase the points behind the Earth, cause this would destroy the topology of the polygons.

Isn't it a simple task? How should I do it?

EDIT - What I've done so far:

The shapefile was downloaded from here to a folder named "shp". The R code is just:

library(rgdal)
countries <- readOGR("shp","TM_WORLD_BORDERS-0.3",encoding="UTF-8",stringsAsFactors=F)
countries <- spTransform(countries,CRS("+proj=ortho +lat_0=-10 +lon_0=-60"))

It gives me the error

non finite transformation detected:
[1] 45.08332 39.76804      Inf      Inf
Erro em .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  failure in Polygons 3 Polygon 1 points 1
Além disso: Mensagens de aviso perdidas:
In .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  108 projected point(s) not finite


texts: Table views from TOC in drop-down list in ArcPy - ArcGIS 10.5.1

texts:
Just upgraded to ArcGIS 10.5.1 and trying to run my self-made Arcpy scripts but there seems to be a problem with the drop-down lists containing table views in this new version. The problem is that a drop-down list icon should appear in the right side of the box that shows all table views loaded in TOC.

However, I've noticed that this icon only appears when there are 1 or 2 dataframes in the document. When there are 3 dataframes (that's what I need and how I've been working with previous versions) the drop-down icon disappears and the only way to go is browsing through folders (folder icon) to find the db view, which by the way slows down a lot the scripts.

This is how it looks like when there are 3 dataframes loaded in the doc:

And this is the right way (seen when removing one data frame from the doc):

I do not include any code because the example itself doesn't require it. It can be reproduced defining a parameter as "Table view" in an empty python script.

I'd like to know if someone else has noticed about this issue and if there is a simple way to fix it.


texts: I have dataset with around thousand geocoded datapoints, each having specific pollutant and its value assigned to it (so it is all in one dataset, some of the datapoints(stations collecting data) are listed several times due to them reporting different pollutants). I would like to check for the spatial correlation of data and literature points me to Moran's I and Getis Ord calculations for this. I believe there are plugins for QGIS that can calculate Moran's I but I did not find one for Getis. Additionally, would somebody be able to explain what is the reasoning for using one or another or is there a need to use both of them when thinking about spatial correlation?


texts: IUniqueValueRenderer on Polyline and Points made my Features all invisible

texts:
I am currently using IUnique value renderer like this:

 foreach (string nodeType in nodes)
      yield return new KeyValuePair<string, ISymbol>(nodeType, CreateNewSymbol() );

public static ISimpleLineSymbol CreateSimpleLineSymbol(IRgbColor rgbColor, Double width, esriSimpleLineStyle style)
{
    SimpleLineSymbolClass simpleLineSymbol = new SimpleLineSymbolClass();

    simpleLineSymbol.Color = rgbColor;
    simpleLineSymbol.Style = style;
    simpleLineSymbol.Width = width;

    return simpleLineSymbol;
}

finally this is called during creation of the layer:

  /// <summary>
        /// Set's the unique symbol types on the feature layer
        /// </summary>
        /// <param name="featureLayer"></param>
        /// <param name="layerType"></param>
        public static void SetFeatureLayerSymbols(IFeatureLayer featureLayer, IEnumerable<KeyValuePair<string, ISymbol>> layerSymbols)
        {
            IFields fields = featureLayer.FeatureClass.Fields;

            //Get indexes of certain fields
            int id = fields.FindField("TYPE");

            IUniqueValueRenderer uniqueValueRenderer = new UniqueValueRenderer();

            //The fields which are used to map to the symbol
            uniqueValueRenderer.FieldCount = 1;
            uniqueValueRenderer.set_Field(0, "TYPE");

            foreach (KeyValuePair<string, ISymbol> values in layerSymbols)
                uniqueValueRenderer.AddValue(values.Key, null, values.Value);

            if (uniqueValueRenderer.ValueCount > 0)
                (featureLayer as IGeoFeatureLayer).Renderer = uniqueValueRenderer as IFeatureRenderer;
        }

I can see the index in the TOC showing the proper unique values and the colors, but on the map it is invisible, I CAN HOWEVER HIGHLIGHT IT.

Thanks and Regards,


texts: QGIS Python Plugin - Prevent Dialog From Closing Immediately after OK Clicked

texts:
I've generated a QGIS plugin using the QGIS Plugin Builder v3.0.3.  One issue I am encountering with it is - as soon as I click the "OK" button (that invokes the run method) - the dialog closes.  Is there anyway to prevent this behavior?  The plugin I am developing is downloading data from a REST api - which can take some time depending on extent.

Ideally I'd like to keep to the dialog open - and showing an indeterminate progress bar - until the data is done downloading.  What happens now is the dialog closes immediately, and then a few seconds later - the raster is added to the map.

Run method from the plugin below.  I didnt alter anything else in the code templates generated by the Plugin Builder.

def run(self):
    # show the dialog
    self.dlg.show()

    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # build dem_url from input form parameters
        # what you see here is just a dummy url for this snippet
        dem_url = 'http://localhost/dem'
        r = requests.get(dem_url,stream=True)
        output_dem = self.dlg.lineEdit.text()
        with open(output_dem,'wb') as out:
             for chunk in r.iter_content(chunk_size=1024):
                 if chunk:  # filter out keep-alive new chunks
                     out.write(chunk)

        iface.addRasterLayer(output_dem, os.path.split(output_dem)[1])


texts: Why does my ITopologicalOperator Intersect result have incorrect area?

texts:
I have used ITopologicalOperator for getting common Area of two polygons but the result is not correct. May anyone please tell me how to get the common area of two features.

//Make Topological operator of geometry here
ITopologicalOperator pTopoOperator = (ITopologicalOperator)pGeometry;
//intersect of the above topo operator with another geometry
IGeometry pGeometry1 = pFeature1.Shape as IGeometry;
IGeometry pGeomResult = pTopoOperator.Intersect(pGeometry1, esriGeometryDimension.esriGeometry2Dimension);
//Typecasted to area viz. common area
IArea pCommonArea = (IArea)pGeomResult;
double dblCommonArea = pCommonArea.Area;

The Area value in dblCommonArea is not coming correct.


texts: Polygon to Raster tool while keeping less than maximum values

texts: I'm trying to convert a very large, detailed polygon to a raster with a coarser resolution in ArcGIS 10.1 using the polygon to raster tool. The value field is a percentage and the cellsize must be the coarse resolution. My problem comes when I'm trying to choose a cell assignment type. If I choose cell-center or no option I get percentages from 73 to 100, while both maximum and maximum combined area give me roughly 50% and higher. I need the raster to include values of 10% and higher, but I don't know how to accomplish this with the options I have. I've tried making the priority field the percentage as well, but it didn't help.


texts: Selecting files with same first 5 characters in filename for merging?

texts:
I have a number of polygon feature classes in a geodatabase dataset. The file names are as below:

e.g. IHO1a_xxx-xxxx_ALB,IHO1b_xxx-xxxx_MBES_10m, IHO1b_xxx-xxxx_MBES_2m

I would like to get a python script which will select all files starting with exactly the same first 14 characters in the filename and merge them.

import os,arcpy

# Define location of .gdb
arcpy.env.workspace = "D:\Users\Documents\H-11_Documents for   Reports_Paul\Test_Model_Builder\ICP_Intermediate.gdb\Deliveryxx_Clipped_Po        lygons"

#Make workspace a variable
workspace = arcpy.env.workspace

List = []

for dirpath, dirnames, filenames in arcpy.da.Walk(workspace, datatype="FeatureClass", type = "Polygon"):
    for filename in filenames:
        if filename.startswith(filename[0:13]):
            List.append(os.path.join(dirpath, filename))
print List
arcpy.Merge_management([filename,filename], 'D:\Users\cadetpn\Documents\H-11_Documents for Reports_Paul\Test_Model_Builder\Temp.gdb\Temp_Data\merge%n%')

##To keep console window open
raw_input("Press enter to exit...")


texts: PostGIS - interfacing with CADCORP 7.1 and QGIS - user 'postgres' can import, but no other roles can

texts:
I'm new to PostGIS (reasonably) - but I'm plagued with one on-going problem:

I'm able to access PostGIS from various GIS systems using the main 'postgres' user.
Using a newly created user 'gis_viewer' - I've granted SELECT and USAGE on SCHEMA etc.
Perhaps the only thing I have not done yet is try to get 'gis_viewer' to INHERIT from a role 'rl_gis_viewer'.

Opening QGIS db manager - I can see that the 'gis_viewer' has access to various schema and rights within those schema.

However when I drag that content into QGIS - I get an error message:

  Layer is not valid -The layer dbname='MasterMap' host=VMSER-ITG-063 port=5432 user='gis_viewer' password='****' sslmode=disable table="garages"."garages_merge" (wkt_geometry) sql= is not a valid layer and can not be added to the map.

The GIS stackexchange also implies that this is a problem with the layer (unique ID issue) -  but that layer is accessible using 'postgres' without these messages. So unless there is a conditional access to layers - I can only assume I've set up gis_viewer incorrectly or PostGIS incorrectly.

Using CADCORP 7.1 - I don't even get that far - apparently the Database does not have PostGIS configured.

This is very much a stumbling block to providing my workplace with access to PostGIS at a SELECT, USAGE level.


texts: I am looking for territories that are divided between two or more countries in a non-contiguous way. E.g, country A receives the western part of the territory, then country B receives a part to its east, then country A again receives a part to its east, and finally country B receives the eastern part. Where can I find such examples?


texts:
I've downloaded the .DEM files from the NASA SRTM30 dataset (found here) however trying to open these yields ERROR 4: "E140N90.DEM" is not recognized as a supported file format.

This is the the same while either using the windows binaries with python bindings and OSGeo4W shell. OSGeo is using GDAL 2.2.2 and USGSDEM files are listed under the supported formats in gdalinfo --formats

 ...
 USGSDEM -raster- (rwv): USGS Optional ASCII DEM (and CDED)
 ...

I've also tried opening the files in QGis, which has yielded the same result. So is there another way to open the files, or convert the files into a supported file format?


texts:
I have a strange behavior with my map.

I wanted to have to basemaps and one overlay.
I've defined classes but when I try to change them via LayerControl I only get one same map.

Could you look at my code?

This script only runs on IE and not on FF or Chromium.

    <!DOCTYPE html>

<meta charset=utf-8 />
<title>A simple map</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

<link href='https://api.tiles.mapbox.com/mapbox.js/v2.1.4/mapbox.css' rel='stylesheet' />


texts:
I have map "my_pointmap" with 3D points with ids=(1,2,3...N). I have also list of subsets of these point ids , S = [(1,2,3,4),(2,4,10,12), ...]. I would love to get GRASS GIS vector map with |S| 3D convex hulls.

for s in S:
    s_ids_string = ",".join(map(str, s))
    gscript.run_command('v.hull',
                     input="my_pointmap",
                     output="my_hullmap",
                     layer = 1,
                     where = "id in ("+s_ids_string +")")

This of course works only for the first iteration. The second one throws error that the layer "my_hullmap" already exists. I could create |S| maps called "my_hullmap_%s". But since S is quite long, I would end up with humongous number of maps.

Note: Since ST_Hull in PostGIS is not doing what it should/could, I am here. Once I will have the map of 3D hulls, I will use v.out.postgis to store it back to postgis.


texts:
The out-of-the-box labeling engine does a fine job with forcing all labels to display despite of overlap. The same cannot be said of Maplex.

I have a polygon layer obtained from the Ontario Provincial government which is solely used for labellings geographical features. (Geographic Name Extent) This is a polygon layer that contains all known "Geographical Name" features along with attributes such as various categories and names in English and French. (Federal government CanVec/NTDB uses a point layer for a similar purpose)

I have a layer that I maintain in which i store the Label schema; all the label 30 definitions and symbols. I simply just import the label schema and turn the labels on. (export and manually correct if needed) I am in the process of updating this layer to work with the Maplex engine but I am unable to place all the labels.

I have 30 label classes
each class is defined by a simple SQL query
each class has a different symbol
each class has different placement properties depending on type
all label classes are set to "Never Remove (Allow Overlap)"
classes with highest priority have a weight of 1000, lowest priority 0

Still many labels are not placed....
Any ideas?

EDIT: I moved the layer to an empty map and I can confirm that when Maplex is off everything is labelled. Turning Maplex on un-places all SMALL POLYGON labels. However, if I zoom extremely close the label appears. I am beginning to feel like Maplex might have been a waste of money. I tried to bring in the layer in question from scratch and creating just one default label class using Maplex - turned on "Never Remove" with same results. All small polygons indicating a location of a town or a bay will not label until over-zoomed. A bug?

EDIT 2: It gets better! I exported all to annotations and out of 1002 features only 429 annotations were created. (Out of which 4 are unplaced) 573 annotations Which I was hoping to "place" manually are missing.


texts: ERROR 6: EPSG PCS/GCS code error when running QgsGeometryAnalyzer().buffer

texts:
I am trying to buffer points from a FileGDB in a PyQGIS script.  I setup my environment with the following batch file and run the python script below.

The problem is that I get ERROR 6: EPSG PCS/GCS code 102184 not found in EPSG support files when QgsGeometryAnalyzer().buffer executes.  The EPSG should actually be 3400 for the input points.  How can I get this to work?

NOTE:
I can enter the commands in the QGIS python console and it works, but the standalone script does not.

setup.bat

SET OSGEO4W_ROOT=C:\OSGeo4W
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
SET QGIS_PREFIX_PATH=%QGIS%

SET PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
SET PYTHONPATH=%QGIS%\python;%PYTHONPATH%
SET GDAL_DRIVER_PATH=%OSGEO4W_ROOT%\bin\gdalplugins
SET GDAL_DATA=%OSGEO4W_ROOT%\share\gdal

buffer.py

import os, sys
from qgis.core import *
from qgis.analysis import QgsGeometryAnalyzer

outputPath = r"D:\buffer"
gdb = os.path.join(outputPath, "data.gdb")

size = 500

QgsApplication.setPrefixPath("C:\\OSGeo4W\\apps\\qgis",True)
QgsApplication.initQgis()
tpts = gdb + "|layername=tpts"
tbuff = outputPath + "\\buff.shp"
tptsQgsLyr = QgsVectorLayer(tpts, "tptsQgsLyr", "ogr")
QgsGeometryAnalyzer().buffer(tptsQgsLyr, tbuff, size,  False, True, -1)
QgsApplication.exitQgis()


texts: I'm working on QGIS, How do I intersect a DEM and a polyline to get the longitudinal section of a river?


texts: "Select a tool" error Saga: Accumulated cost (isotropic)

texts:
I'm trying to run Saga's accumulated cost operation in QGIS on a DEM and set of raster destination points, but keep receiving an "Error: select a tool" in the log after a failed run. There are no other options to specify in the Saga Accumulated Cost algorithm (at least that are presented in it's processing toolbox dialogue), and am not sure what is producing this error. How do I resolve this?

Alternatively, is there another means of generating an accumulated cost map? Saga's algorithm is proving to be extremely difficult to use.


texts: IFeatureBuffer - performance considerations (ArcObjects)

texts:
I'm trying to optimize my code to make it faster. I have a loop that reads some string input and creates points with the help of IFeatureBuffer and IFeatureCursor.

If I create the buffer inside tho loop, it slows the code down.

For i = 0 To Counter
  pFBuffer = pFClass.CreateFeatureBuffer
  ...
  pFBuffer.Value(iField) = ...
  pFCursor.InsertFeature(pFBuffer)
Next

If I move pFBuffer = pFClass.CreateFeatureBuffer above the loop (i.e. create the buffer once only) it works significantly faster. However, I have to "clear" the buffer at the end of each iteration, as in some cases it stores values from previous iterations.

So the question is: how can I clear the buffer, so that I won't need to create it at each iteration? Are there any other drawbacks or things to consider?


texts: Clipping vector shapefiles more efficiently using QGIS?

texts:
I have a river shapefile, got from r.watershed module, and it looks like this：

I have eight points on the river which I want to analyse, and I only want this river segments：(the second photo below is the river I want)

So I open its attribute table to delete segments I don't want, but I think it's not efficient.

I have tried the intersect tool in Geoprocessing tools, but the results were not what I want.

Is there more efficient way to achieve this goal？


texts:
I have high res raster data covering 2million km sq, with a field where each cell has a value from 0 to 100. I need to create a "mask"  using the cells with a value over 12 in that field. I then need to clip vector data to that mask.
I was wondering if there is a good method for selecting cells from rasters and using that selection (with reduced noise and smoothed edges) to clip vector data?
I don't have spatial analyst extension. I am using ArcGIS for Desktop 10.1.

-


texts: Calling help for ArcPy function from within IDE

texts:
In IDLE (for Python v2.7.2 within my ESRI path), when I type:
help(arcpy.gp.Slope_sa)

I get:
Help on function <lambda> in module arcpy.geoprocessing._base:

<lambda> lambda *args

which is not particularly helpful.
How do I access proper help?


texts: Iterating through geodatabase with multiple feature datasets using ArcPy?

texts:
I am trying to create a Python script that batch projects every shapefile in a geodatabase to a temporary geodatabase that I created. I realized that this script that I created does not iterate through the multiple feature datasets of the geodatabase.

The code that I created only projects the file "breaklines". I need to also iterate through "Contours" and "Hydro". I know that I need to create a double loop but I am not sure how to do it. Here is the code that I started below

import glob, os, shutil, time, arcpy

path = r"C:\...\SARA_New_Contours.gdb"
temp = r"C:\...\Temp.gdb"

string1 = "\\metersproj-"

arcpy.env.workspace = path
out_coor_system = arcpy.SpatialReference("C:\...\NAD 1983 (2011) StatePlane Texas S Central FIPS 4204 (Meters).prj")

for fc in arcpy.ListFeatureClasses():
    prefix = string1 + fc
    out_dataset = temp + "\\" + prefix
    arcpy.Project_management (fc, out_dataset, out_coor_system)

A side question: How would I just loop through one dataset, for example: just the feature dataset called "Contours"?


texts:
When I go to Raster -> Extraction -> Clip Raster by extent, I get the error "/bin/sh: gdal_translate: command not found ".  I am using QGIS version 3.0 on a MacBook Pro.

Anyone know how I can fix this?


texts: Populating SQL options for parameter in ArcGIS tool?

texts: In ArcGIS Desktop 10 I am trying to create a tool where individuals can select a variety of records using a SQL expression that is then passed to a calculation. I can add a parameter under the tool, and select SQL expression, but when I click on the SQL expression box its not "populated" by potential fields, and therefore unable to get unique values. I thought this might be because the tool needed to know which .shp file it was calling, so I added a higher level parameter to specify the .shp file, but still no luck, the SQL parameter remains entry.


texts:
Link

This should be the Topographic Position Index

Pls run the Code and Zoom in and out

you will notice that the representation changes for different zoom levels, also Values from the Inspector Change

seeking for clearence


texts: Detecting CTRL key in OpenLayers select interaction?

texts:
I have an OpenLayers select interaction with an add condition event handler.

const addCondition = (event: MapBrowserEvent):boolean => {

  return true;
}

In this handler I need to determine if the CTRL key was pressed when the user clicked on a map feature. It does not appear the MapBrowserEvent contains the information to determine if the CTRL key was pressed. Is there a way to do this?


texts: QGIS as a pre and post processor for modflow

texts: I want to find out how QGIS can be used as a pre-processor for modflow, i.e delineating the model boundary and boundary conditions etc. exporting to modflow for analysis, viewing and/or mapping the results in QGIS. I am new to QGIS and GIS in general but I was told that this is possible in QGIS.


texts: Calculating density using information linked to polygons?

texts:
I'm fairly new at GIS and am using ArcMap to calculate the population density of lizards on an island. On our island, we have about 80 trees, and each tree is mapped as a polygon on a single layer. We have a database of several thousand lizards that we have captured on these trees (and noted capture location). Our goal is to add the lizards to the map to calculate density across the island. Right now, we have attributes for each tree that are #males, #females, #total, which we hoped would give use the ability to calculate both density and sex ratio.

Can we use these attributes to calculate density of specific areas with a given kernel size? Based on what I've read, it seems like you might only be able to calculate density of lines and points.


texts: How to automate drawing a polyline from 2 intersecting datasets based on proximity?

texts:
I've been stuck on a situation recently where my company's data has to be integrated with another.
Basically there are two point shapefiles comprising of >50000 points, and as shown in the image below, one data set lies between the other spatially.
The output I need would be around <100 polylines based on a unique identifier within the two datasets.
The problem occurs because once they are merged, the FID number (which I would usually use on the Line type within the points to line tool) no longer has any spatial use as my own dataset is at the top of the attribute table and the other at the bottom (no longer in spatial order).

Can anyone think of a solution to this?
A polyline tool based on proximity would be perfect for this situation but I am unaware of one existing...
(I've already thought about sorting by northing/easting, but that creates problems when the data does not go in a continuous orientation)


texts: Making legends for each layers and switching between them using Mapbox GL JS?

texts:
I created a map using Mapbox GL JS with several layers using this example

At the same time I don't have any spatial knowledge in Javascript. It's restrict me to create legends for each layers. The particular legend(s) must show when particular layer(s) switched on.

I find function L.mapbox.legendControl but how can I use it?


texts: Geoserver: huge CPU load when delivering tiles at higher zoom levels

texts:
We are working with geoserver for the first time. We have it running on a quite powerful server, 32GB and 12 kernels.

For the raster we want to deliver we have already seeded all cache tiles for all zoom levels, i.e. from 0 to 19, which is the max zoom.

When looking at the CPU load while accessing the tiles there is no unusual behavior of the CPU load (around 10% to 30%).

However, when one gets to zoom level 19 the CPU load jumps to 300+%, reaching even >1000% (corresponding to the utilization of almost all available kernels).

How can this be? There is no tile generation, since they are already all available, so the cpu should not be used in this manner. Their delivery should not be different from other zoom levels, as far as I understand it, so there should not be more CPU activity than with lower zoom levels.

Is this a bug?

Does someone know the reason and maybe also the solution to this problem?

I prepared a video to show what a mean, please look at https://www.youtube.com/watch?v=f0VbULDyj7c


texts:
Trying to use raster extraction clipper to prepare a selected area of the map for generating contours - keep getting the error message:

  ERROR 1: TIFFOpen:everest_gmted30.tif: Permission denied

Any clues on what to do about this?


texts: Get the bounding box for each feature in the source layer in QGIS

texts:
I have a polygon layer with multiple polygon features drawn on it.
I need a tool that will draw a bounding box around each polygon feature individually. Right now everything I try in the tool box draws a single bounding box that encompasses all the polygon features. Is there such a tool?


texts: Icons in ArcToolbox window acting very strange?

texts: I am running ArcGIS 10.3.1 for Desktop and for some reason my icons in the ArcToolbox and in other windows are behaving very very strange.  Sometimes no icon at all appears, other times the completely wrong icon appears. The icons will randomly change each time I load a new MXD.  I have searched the net for an answer but nothing comes up.  I attached a screenshot below.  First off anyone know why this is happening and if this is a serious problem that would require a re-install or other corrective measures?


texts: How to add a simple WFS layer from GeoServer to OpenLayers map?

texts:
I'm trying to add one of the included datasets with GeoServer to an OpenLayers map. I've followed the OpenGeo tutorials, seen similar thread on here, but I can't seem to crack it. Can somebody take a glance at my code and GeoServer setup and tell me where I'm going wrong?

Here's my code:

//WMS map
world = new OpenLayers.Layer.WMS("Global Imagery", "http://maps.opengeo.org/geowebcache/service/wms", {
    layers : "openstreetmap",
    format : "image/png"
});
map.addLayer(world);

//WFS
parks = new OpenLayers.Layer.Vector("WFS", {
    strategies : [new OpenLayers.Strategy.BBOX()],
    protocol : new OpenLayers.Protocol.WFS({
        url : "http://localhost:8081/geoserver/wfs",
        featureType : "medford:parks",
        featureNS : "http://medford.opengeo.org/medford"
    })
});

map.addLayer(parks);

map.zoomToMaxExtent();

GeoServer details:
workspace name = medford, Namespace URI = http://medford.opengeo.org
For the feature type I have tried "parks" and "medford:parks"

The layer is enabled, advertised and has a 30 feature return limit.

I'm struggling to put a simple map together.

Updated Solution. Working Now.

wfs = new OpenLayers.Layer.Vector("Fields_WFS", {
        strategies : [new OpenLayers.Strategy.Fixed()],
        protocol : new OpenLayers.Protocol.WFS({
        version : "1.1.0",
        url : GEOSERVER_HOST + ":" + GEOSERVER_PORT + "/geoserver/wfs",
        featurePrefix : "rpid",
        featureType : "FIELDS_SUBSET_BNG_OSGB36_1",
        featureNS : "<namespace>",
        //geometryName : "GEOM", type "Geometry"
        srsName : "EPSG:27700"
}),
    renderers : renderer
})

Does not need to specify the geometry Name in Layer definition.
Since Native SRS was null in our dataset, if we keep geometryName than it takes default, but without it, it uses map projection(27700- in our case).


texts:
On aviation maps (ICAO map 1:500.000) restricted areas are indicated by a line with hachure pointing to the inside of the polygon.

Refer to ICAO Annex 4 Appendix 2 page APP 2-9:

https://www.casa.gov.au/sites/g/files/net351/f/_assets/main/lib91135/an04_cons.pdf

Is it possible to create a line or fill style with QGIS which represents this hachure?


texts: Searching features by text and Zoom to extent of that feature in Cesium

texts: I am using Cesium and I am new to it. I want to implement searching by certain attribute for example Name of the point. The result obtained should be the map zoomed into the point found. How can this be implemented.


texts: solve Topology error in Grass using QGIS?

texts:
I want to build a raster from a vector layer by using "v.to.raster.attributes" in QGIS but its not working and give me an error to build a topology using v.build. I build it but after that i have no layer to add in qgis using grass plugin. it doesn't give me any error. I am not getting what's actually problem is ?

Out put after adding vector layer
"v.in.ogr dsn=D:/Thesis/Data/Mangla_Study_Area/mangla_area_p.shp output=Mangla snap=-1 min_area=0.0001 -o

Over-riding projection check

Layer: mangla_area_p

Counting polygons for 1 features...

Importing map 1 features...

-----------------------------------------------------

Building topology for vector map ...

Registering primitives...

1 primitives registered

54951 vertices registered

Number of nodes: 1

Number of primitives: 1

Number of points: 0

Number of lines: 0

Number of boundaries: 1

Number of centroids: 0

Number of areas: -

Number of isles: -

-----------------------------------------------------

￼Cleaning polygons, result is not guaranteed!

-----------------------------------------------------

Break polygons:

-----------------------------------------------------

Remove duplicates:

-----------------------------------------------------

Break boundaries:

-----------------------------------------------------

Remove duplicates:

-----------------------------------------------------

Clean boundaries at nodes:

-----------------------------------------------------

Merge boundaries:

-----------------------------------------------------

Change dangles to lines:

-----------------------------------------------------

Remove bridges:

-----------------------------------------------------

Building topology for vector map ...

Building areas...

1 areas built

1 isles built

Attaching islands...

Number of nodes: 1

Number of primitives: 1

Number of points: 0

Number of lines: 0

Number of boundaries: 1

Number of centroids: 0

Number of areas: 1

Number of isles: 1

Number of areas without centroid: 1

-----------------------------------------------------

Find centroids for layer: mangla_area_p

-----------------------------------------------------

Write centroids:

-----------------------------------------------------

1 input polygons

Total area: 2.26321E+008 (1 areas)

Overlapping area: 0 (0 areas)

Area without category: 0 (0 areas)

-----------------------------------------------------

Building topology for vector map ...

Registering primitives...

2 primitives registered

54952 vertices registered

Building areas...

1 areas built

1 isles built

Attaching islands...

Attaching centroids...

Number of nodes: 2

Number of primitives: 2

Number of points: 0

Number of lines: 0

Number of boundaries: 1

Number of centroids: 1

Number of areas: 1

Number of isles: 1
Successfully finished"

and after applying topology V.build

"v.build map=Mangla@Mangla_reservior option=build error=Mangla_topology

Building topology for vector map ...

Registering primitives...

2 primitives registered

54952 vertices registered

Building areas...

1 areas built

1 isles built

Attaching islands...

Attaching centroids...

Number of nodes: 2

Number of primitives: 2

Number of points: 0

Number of lines: 0

Number of boundaries: 1

Number of centroids: 1

Number of areas: 1

Number of isles: 1

Building topology for vector map ...

Registering primitives...

0 primitives registered

0 vertices registered

Building areas...

0 areas built

0 isles built

Attaching islands...

Attaching centroids...

Number of nodes: 0

Number of primitives: 0

Number of points: 0

Number of lines: 0

Number of boundaries: 0

Number of centroids: 0

Number of areas: 0

Number of isles: 0
Successfully finished"


texts: Dissolve with common value in a field and not common in another

texts:
I have a shapefile in QGIS and I want to dissolve it with the options:

A field named 'specialid' has the same value

A field named 'city' has different value

Ok the first is very easy. It is basically what dissolve is meant to be
So all same specialid will be unified. But I want those with the same 'specialid' and different 'city' value NOT to be dissolved.
How can I process the second option? Is there a way to select it or do it another way?


texts: Sort PDF Files based on last 7 Characters

texts:
I have a Python script that creates PDF Map Books from pre-existing map PDFs, created in ArcGIS Pro. The PDF file names end with the month abbreviation and year they are associated with (example: CACumulativeSumWithin3000_SanBruno_Dec2017.pdf). I want to order the PDFs in the Map Book based on these final 7 charcacters, but the problem is it orders the PDFs alphabetically, so Apr2018 is first instead of Jan2018.

How can I order these PDFs using the month abbreviation and year?

Here is my code:

import arcpy
import os

arcpy.env.workspace = r"E:\arcGIS_Shared\Python\Export\Cumulative_CA2"
path =  r"E:\arcGIS_Shared\Python\Export\Cumulative_CA2"

### Create PDF Map Book File
CApdfPath = os.path.join(path, 'SanBruno_CumulativeHeatMaps.pdf')

CApdfDoc = arcpy.mp.PDFDocumentCreate(CApdfPath)

### List all files in the output folder
CApdfList = arcpy.ListFiles( "CACumulativeSumWithin3000_SanBruno_*")

### Add market title page to the newly created PDF document
CApdfDoc.appendPages(os.path.join(path, "CATitlePage.pdf"))

### Append market specific pages to the PDF
for CApdf in CApdfList:

    print(CApdf)

    CApdfDoc.appendPages(r"E:\arcGIS_Shared\Python\Export\Cumulative_CA2\%s"%CApdf)
    print(CApdf)

CApdfDoc.saveAndClose()


texts: Losing pixels values when counting by value

texts:
I have a script that calculates number of pixels with NDVI values that are higher than 0.3, 0.5 and 0.7.

My code runs with no errors, but it get value 0 and doesn't really count.

what can be the reason?

var meanImage=ee.Image(byYear.get[0]);

var gt03 = meanImage.gt(0.3).selfMask().rename('NDVI_gt03');
var gt04 = meanImage.gt(0.5).selfMask().rename('NDVI_gt05');
var gt05 = meanImage.gt(0.7).selfMask().rename('NDVI_gt07');

Map.addLayer(gt03,NDVIcolor,'gt03');
Map.addLayer(gt04,NDVIcolor,'gt05');
Map.addLayer(gt05,NDVIcolor,'gt07');

//count the number of total pixels
var c03 = gt03.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: geometry,
  crs: 'EPSG:4326',
  scale: 30,
  });

var c05 = gt05.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry:geometry,
  crs: 'EPSG:4326',
  scale: 30,
  });

  var c07 = gt07.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry:geometry,
  crs: 'EPSG:4326',
  scale: 30,
  });
// print the dictionary output
print(c03,'Total Pixels gt 0.3 ');
print(c05,'Total Pixels gt 0.5 ');
print(c07,'Total Pixels gt 0.7 ');

I have tried also to use listOfImages but still got 0.

My end goal is to calculte the number of pixels with the given values for one image.

link for my code:
https://code.earthengine.google.com/559d0922b50d67cf01e4d5b03a21666d


texts:
I know there already are posts about this subject, I tried to derive a solution for my own problem but the code for field mapping is very complex to me.

I have a Raster 100x100m and point layers that I want to join using a python script. I need the SUM of a field called "weight" (point layer) added to the grid, but I don't know how to define the field mapping.

So far the script looks like this (but doesn't sum up my weight field):

import arcpy
from arcpy import env

env.workspace = r"C:\Users\Daniela\Desktop\test"

fc_list = arcpy.ListFeatureClasses()

for fc in fc_list:
    arcpy.SpatialJoin_analysis("grid100.shp",fc,"C:/Users/Daniela/Desktop/test/test" + "_" + str(fc),"JOIN_ONE_TO_ONE","KEEP_ALL","sum","INTERSECT","","")


texts: Using SelectLayerByLocation and GetCount to Populate a new Field using Python

texts:
I'm attempting to create a script which assesses the intersect between a layer "BldgFootprints" and several polygons "FireBoxMap_1"-"FireBoxMap_44". And from this intersects performs a getCount on a particular attribute of the BldgFootprint layer and uses it to populate a newly created field in the FireBoxMap polygon attribute table. Essentially I'm trying to create a count of single family homes within particular fire districts using a buildings database and then adding this count to the fire district attribute table to tell at a glance how many homes are in each district. I am using PyCharm IDE and ArcGIS Pro.

So far I have:

import arcpy
try:
    arcpy.env.workspace = (r'L:\ProgrammingPro\FireDepartment.gdb')
    fcs = arcpy.ListFeatureClasses("FireBoxMap*")
    print (fcs)
    lyr = arcpy.ListFeatureClasses("B*")
    print (lyr)
    for fc in fcs:
        inFeatures = "fc"
        fieldName1 = "Single"
        fieldName2 = "Multi"
        arcpy.AddField_management(inFeatures, fieldName1, "SHORT", "","","","","","","")
        arcpy.AddField_management(inFeatures,fieldName2, "SHORT", "","","","","","","")
    bldg_lyr = arcpy.MakeFeatureLayer_management(b*,r'in_memory\bldg_lyr')
    firebox_lyr = arcpy.MakeFeatureLayer_management(FireBoxMap*,r'in_memory\firebox_lyr')
    firebox_bldg = arcpy.SelectLayerByLocation_management(bldg_lyr,'HAVE_THEIR_CENTER_IN',firebox_lyr*)

I have a failed to execute error in the MakeFeatureLayer section and don't know where to go from there.


texts:
I'm trying to clip a large polygon shapefile (nearly 440,000 features, 0.5GB) to a shapefile containing a single polygon. Both shapefiles use the same CRS (EPSG:27700).

With QGIS, using Vector > Geoprocessing Tools > Clip, I left it running for nearly 24 hours and it was still at only 7% complete.

I tried the same clip using ArcMap and it finished in 4 minutes.

I tried it again in QGIS, and also tried GDAL/OGR > Geoprocessing > Clip vectors to polygons, but it's been running for an hour and is clearly going to take an eternity.

I'm new to QGIS and am wondering if this typical of its performance when clipping large vector datasets? Or is there a trick to get it perform somewhere close to ArcMap's speed?


texts: Does a high number of SDE_state_lineages slow performance?

texts:
In our SDE database, the SDE_state_lineages table has around 200,000 records, the SDE_states has around 900, and our SDE_versions table has around 1,300 records.

We've been having performance issues lately. The layers are very slow to draw. We compress daily, but we aren't able to reconcile the existing versions.

Is this high number of records in the SDE_state_lineages causing poor performance?


texts: Place a special function or action with QGIS

texts:
I use QGIS and  want to place an option that when pressed, I perform a filtering without the need to enter commands.

I am currently using the program to carry out an electricity project, where a series of high and low voltage electric poles and networks is loaded on a map, but when I finish loading these poles electric and networks layer I must correct the mistakes made in the tracing of networks or filling of spreadsheets, as we are several operators, at the end of the afternoon you must use the dialog expression to filter by operator the mistakes you make and I must place the command again and again, so I would like to add a button or option that allows me to select an already established filter.

How could this be done?


texts:
I am a fairly strong user of ArcGIS, having started with version 8 long ago. But I'm looking for an inexpensive alternative right now for a client, and I want to investigate QGIS as a possible solution.

Thing is, I have never touched it before. I see mention of Grass databases and OSGEO4W shells and Python. I don't know what any of that means. No clue.

So where can I find a beginner's guide to QGIS, that more or less walks a newbie through the various components, how they work, why they're needed, and gets me to the point that I can import a .gdb and make some sample maps and queries to show my client?


texts: Is it possible to perform Interpolation in QGIS with input as a vector layer which has a set of points collected following a specified path?


texts: Convert List by Source to List by drawing order using arcpy

texts: When I do some calculations, it results in adding a table to TOC which converts "List by drawing order" to "List by Source". I want to convert it back to "List by drawing order" using arcpy and I am not able to find anything relevant to it in help files. Any suggestion?


texts: Data-Driven Styling (Mapbox GL JS) for Choropleth Maps - Using Multiple Fields

texts:
I have been diving into the Mapbox GL JS library lately and am super impressed and excited by the possibilities and capabilities. Currently I am looking for a solution to building choropleth maps, each based on a different indicator field which could be selected by toggle-button or dropdown.

I've come across an example that is similar to what I have in mind, but I am wondering if this might be done more efficiently with the advent of data-driven styling.

Another goal would be to generate natural breaks values (or other classification method) dynamically when an indicator field is selected based on the range of data, and then push these values to a legend.

Finally, joining data to a google sheet for field/indicator values would be amazing, but I don't know if this sort of join of the data to a google sheet is possible for polygon data. However it is done with X, Y point data in several examples here using util.getJSON.

Data-driven styling for fills example.


texts: Styled markers : how to get bouncing markers with location dots?

texts: I want to mark important location with bigger markers (pointer with forks) and smaller location with normal markers (colored dots). I want big markers to bounce on hover. How to achieve this effect in Google Maps or using Leaflet.


texts: GDAL VRT Format - do I understand the concept?

texts:
Okay, so I'm a little fed up with different raster-in-database solutions that I tried over the last 2 weeks which don't work currently (PostGIS raster) or can't handle my volume of data (Rasterlite), so I'm reading about the VRT format. Bear with me, I'd like your opinions on this, and also if what I propose makes sense at all:

VRT is basically a metadata XML file describing various properties of the actual raster file, like pixel dimensions, geolocation, etc..

Step 1: I converted (not really, just generated VRTs) all my GeoTIFF files into VRT.

Now, since I see I can merge the VRT, I assume its possible to create a mosaic of VRT metadata files, one big VRT to hold all metadata. Great, that means, metadata will probably be used for optimizing the rendering since QGIS will decide (assumption) which tiles fall into the current canvas extent and which not...

Step 2: create a VRT mosaic

Okay, now the next benefit I see is, I only need to add one entry into the QGIS table of contents and I can set different things upon it, like scale dependant rendering and null data value transparency if needed. One could argue, you can simply load GeoTIFFs into QGIS and put them in a group, true, but you don't have any control over a group except for setting the common SRS.

If I wanted to see my orthophoto on a small scale (large zoom-out) and still have speedy rendering I would repeat the steps, only with a resampled set of original raster and end up with a few big VRT datasets.

Does this make any sense? :)


texts: Drawing arrow on line using SLD of GeoServer?

texts:
I have a line layer in geoserver as in the 1st image.

How can I draw an arrow on the line like in the 2nd image?

The line.rar includes a shapefile and a sld file.

Can edit a sample for me?

line.rar(download)


texts: convert list of points to line or polygon

texts:
I have a list of Points in a PostGIS table. The rows can be grouped by a date

something like:

myTable
 id INT
 date DATE
 lat NUMBER
 lng NUMBER
 geom GEOMETRY(POINT)

now I need a table where all points of one day are collected in a new geometry (sorted by id)

targetTable
 id2 INT
 date DATE
 geom2 GEOMETRY(POLYGON)

how can I do this with PostGIS or GeoKettle, I cannot use any other client since this transformation has to be done in batch?


texts: Making fields conditional using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop to build a large inventory of urban trees to be used for tree risk assessment purposes. This inventory was originally in a relational Microsoft Access database, and I have spent the last few months structuring a new database using ArcGIS. I have created relationship classes between the main Feature Class ("Tree_Info") and other secondary tables, such as the "Inspection" table with risk assessment information.

In the "Inspection" table, there is a field called "Significant" which is a boolean 1-Yes 0-No. Not surprisingly, the "Significant" field represents trees of significant importance.

I have 4 other fields in the "Inspection" table which I ONLY want populated when "Significant" is populated by 1-Yes. These fields include "Rem_Mod" (Can the tree be removed or modified? Yes or No), "Sig_Type" (Signficance Type: historic, cultural, memorial, etc.), "DBH" (Diameter at Breast Height), Crown_Sprd (Crown Spread). If the "Significant" field is populated by 0-No, then these 4 fields can be left null.

I've looked into Subtypes, and I'm not sure this is what it is intended to do.

I'm not a programmer, and I am relying on the native ArcGIS tools, unless someone has a different simple solution. I'm at a loss of what to do. Does anyone have any suggestions on how to make these 4 fields conditional upon the tree's significance?


texts: Why does the difference between two overlapping rasters exist?

texts:
I have two similar rasters of the Mount Vesuvius: one spans 20 kilometers, and the other one 50 kilometers. Both rasters are in WGS84 geographic coordinate system, and they overlap between each other:

When I use QGIS Raster calculator to calculate the difference between them, a new raster file is created which is totally black - so it shows that there is absolutely no difference in pixel values on the part where they overlap:

However, when I project both rasters to Azimuthal equidistant projection, and again check for the difference with Raster calculator, the difference between these two reprojected rasters exists! Here is how it looks like:

Why is this so?

The difference does not exist between the two rasters when they are in WGS84 geographic coordinate system, but once they are projected to Azimuthal equidistant projection, a significant difference appears?

Why?

I used the cubic resampling method, but I got the similar upper difference with the bilinear one.

Just for the sake of checking, I tried using Transverse Mercator projection (instead of Azimuthal equidistant), and again a difference exists between initial two projected rasters. Why is that so?

Additional information:
To project the initial WGS84 raster files, I used the following syntax:

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_ts=40.81266 +lon_0=14.414252" -r cubic -of GTiff "C:/vesuvius_radius_20KM.tif" "C:/vesuvius_radius_20KM_cubic_aeqd.tif"

To project them to Transverse Mercator:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32633 -r cubic -of GTiff "C:/vesuvius_radius_50KM.tif" "C:/vesuvius_radius_50KM_cubic_tm.tif"

Here are the initial raster files (in WGS84):

https://www.dropbox.com/s/98hcftjsnmrqs8p/vesuvius_radius_20KM.tif?dl=0 https://www.dropbox.com/s/a68921e6tpszt0f/vesuvius_radius_50KM.tif?dl=0

Here are their Azimuthal equidistant projected rasters:

https://www.dropbox.com/s/4mexe0rmieam9ri/vesuvius_radius_20KM_cubic_aeqd.tif?dl=0
https://www.dropbox.com/s/0hcle0b151kbod1/vesuvius_radius_50KM_cubic_aeqd.tif?dl=0

And the difference between two Azimuthal equidistant projected rasters:

https://www.dropbox.com/s/2w39hgb9kcrzb8h/vesuvius_radius_20-50KM_cubic_aeqd.tif?dl=0


texts: Reduce number of vector field arrows in qgis

texts: OK, now that I have figured out how to get my vector fields in qgis, I can't figure out how to reduce the number of vector arrows on the map. Is there a setting for how many vectors are sampled? Thanks.


texts: I ma currently developing a web mapping app. How do I make it possible for my app to clip and save area of interests or area I am working on. I can create a vector layer, modify it and delete if I want to, but how can I clip using OpenLayers 2.


texts: Map Service Works in Flex, Won't Display using JavaScript

texts:
As the title states, a map service published on my server displays and can be identified just fine when added to a web map in Viewer for Flex. I am slowly learning to use the API for JavaScript. My current JavaScript code works as intended when using a third-party map service for testing purposes. When I replace it with the service I published, it does not display the service. The rest of the map displays correctly, just not my service. I cannot figure out why it works in Flex, but not JavaScript. Any one run into this or have any ideas?

Thanks!

EDIT: Map service, dynamic. Arc 10.1, JS API 3.10

The service does display in the above example, and when you view it through the REST Services directory.

The error I'm getting is "Unable to draw graphic (geometry:null, symbol:null): null" ... but I am not finding any null geometry.


texts: Import shapefile is not working in pyshp?

texts:
I am just trying to use "pyshp" package to open shapefile in Python and use the Attribute data for further development.

Firstly I have downloaded pyshp package for python 2.7 and installed it in my computer.
But, whenever I want to import shapefile, I got following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: No module named shapefile


texts:
I am trying to get the nearest node from a point(Id 1) . The coordinates (lat/long) are input by an user . I am using some functions like pgr_nearest_node() or get_nearest_node() from get_nearest_node function reference

1: http://www.sqlexamples.info/SPAT/postgis_nearest_point.htm. The Problem Is that the point take the nearest node using straight line between them(928-3670-317-929 Id nodes). How can i get the nearest node in the network ? . I am thinking about to use pgr_dijkstra function but it need to a source id so i return to functions commented before.


texts: What does number placed next to bar chart in Table of Contents by ArcGIS for Desktop mean?

texts:
Working with bar charts in ArcGIS 10.2.1

What does the 0.22 represent on the chart symbology in the Table of Contents?


texts: sf equivalent of R maptools package's SpatialLinesMidPoints

texts:
I've been moving most of my workflow to using the sf package in R, and I'm trying to find the equivalent to the maptools package SpatialLinesMidPoints function - I've used this function frequently over the years and had modified to get line endpoints as well (here).

Is there an equivalent function using sf package?  I've looked at st_cast and st_line_sample, but want to explicitly get midpoints of lines. I'm guessing perhaps I can make use of either st_cast or st_line_sample with slight modification to get just midpoints but struggling a bit with approach.

Just to be explicit, I'm looking for equivalent of doing:

library(spData)
library(maptools)
library(sf)

seine <- as(seine, 'Spatial')
seine_midpoints <- SpatialLinesMidPoints(seine)
# To verify:
plot(seine, axes=T)
plot(seine_midpoints, add=T, col='red')


texts: Does File Geodatabase performance degrade as it fills up?

texts:
I'm having a memory leak issue. Upon investigating this more, it seems the file based geodatabase I write to in a loop, is growing very large, and as it grows, it significantly degrades the performance of the scripts I am running.

Any ideas how to optimize the configuration of the fgdb?  Or how to speed it all up?  I am not writing 'in_memory', I am using aggregatePoints to create a temp featureclass (which I delete), and I buffer this FC, which I keep.

However, it seems to get slower, and slower, and slower...

def createGeom1(geom, scratchDB):
    filetime = (str(time.time())).split(".")
    outfile = "fc" + filetime[0] + filetime[1]
    outpath = scratchDB + "tmp1.gdb/Polygon/"
    outFeatureAggClass = outpath+outfile + "_Agg"
    arcpy.AggregatePoints_cartography(geom, outFeatureAggClass, "124000 meters")

geom is a collection of Points, scratch the scratch area (local gdb I am using).

Just looping through a list of files, I call a procedure that creates  alist of geoms (and doesn't degrade) and then call this.  Doing this, will see this function, createGeom, degrade significantly, and the previous one, not a bit.


texts: How to reproject a raster to Winkel Tripel projection?

texts:
I'm trying to project a raster map (in ASCII grid format) from longlat to the Winkel Tripel projection. But I'm running into errors, for instance gdalwarp won't accept the proj-string.

$ gdalwarp -t_srs '+proj=wintri +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0' popd.asc popd_wintri.tif
ERROR 1: Translating source or target SRS failed:
+proj=wintri +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0

And if I try it in GRASS I get the following output

GRASS 6.5.svn (Wintri):~ > r.proj input=popd_2000AD location=Latlong output=test22 --o
Input Projection Parameters:  +proj=longlat +no_defs +a=6378137 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0
Input Unit Factor: 1
Output Projection Parameters:  +proj=wintri +lat_0=0.0000000000 +lat_1=0.0000000000 +lon_0=20.0000000000 +a=6378137 +rf=298.257223563 +no_defs
Output Unit Factor: 1
WARNING: pj_transform() failed: non-convergent inverse meridional dist
(...)
ERROR: Input raster map is outside current region

I got the same error when I tried projectRaster() in R:

Error in projectRaster(test, wintri) :
  error in pj_transform: non-convergent inverse meridional dist

On the GRASS website there's a screenshot of a raster map in Winkel Tripel-projection, so I'm guess it's possible to accomplish this in GRASS. Any help would be appreciated!


texts: PGrouting encapsulated queries error [quotes problem]

texts:
I'm facing a quote problem on a pgrouting query.

What I'm thinking is wrong is that my query only works if I use one, then four, then eight quotes on encapsulated queries.

Is there a way to simplify this ?

SELECT *
FROM pgr_pointsaspolygon(
'SELECT id::int4, st_x(ways_vertices_pgr.the_geom) as x, st_y(ways_vertices_pgr.the_geom) as y
  FROM pgr_drivingDistance(
    ''''SELECT gid as id, source, target, ST_Length_Spheroid(the_geom, ''''''''SPHEROID["WGS 84",6378137,298.257223563]'''''''') as cost FROM ways'''',
    7,
    550,
    false,
    false)
INNER JOIN ways_vertices_pgr on (ways_vertices_pgr.id = id2)');

Thanks in advance


texts: Converting coordinates into Longitude/Latatitude form that CartoDB understands?

texts:
Does anyone know how to convert this infomation into longitude/latitude coordinates for georeferencing on cartodb?

UTC: Posn: 03:00N – 105:10E
LT: Posn: 01:28.0S – 116:48.0E
LT: Posn: 08:44.7S - 013:17.9E

These data represent locations in the middle of oceans.

Is there a simple rule/website for conversion?


texts: Summing attributes of different shapes into another (height of contours + height of buildings)

texts:
I'm using QGIS 2.18.

I need to calculate the height of the rooftop of each building in a neighborhood. I have a vector shape of the buildings, with the building height attribute (from the ground to the top) and the shape of the contours (the height of the ground itself).

I need the information of the height of the ground of the centroid of each building and sum it into the height of the building in order to have the total height of the rooftop, in a shape of the buildings.

Here's a link to what I'm trying to do.


texts:
I'm using emissions data in the form of netCDF files from the CMIP5 project.  I've loaded this data into R as a rasterBrick and can see that it has 6 layers.  My goal is to be able to use this data with Google Earth maps.

My understanding is that plotKML does this.  However, I've had a difficult time getting this to work smoothly.  After contorting paths, I had some success:

> plotKML(bdat[[1]])
Plotting the first variable on the list
KML file opened for writing...
Writing to KML...
Closing  bdat__1__.kml

This produced a file:

<kml xmlns:xsd="http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd"
xmlns:xmlns="http://www.opengis.net/kml/2.2/" version="1.0">
  <Document>
    <name>bdat__1__</name>
    <open>1</open>
    <Folder>
      <name>SpatialGridDataFrame</name>
      <GroundOverlay>
        <name>colour</name>
        <altitude>0</altitude>
        <altitudeMode>clampToGround</altitudeMode>
        <Icon>
          <href>X54786.png</href>
        </Icon>
        <LatLonBox>
          <north>90</north>
          <south>-90</south>
          <east>180</east>
          <west>-180</west>
        </LatLonBox>
      </GroundOverlay>
    </Folder>
    <ScreenOverlay>
      <name>Legend</name>
      <Icon>
        <href>X54786_legend.png</href>
      </Icon>
      <overlayXY x="0" y="1" xunits="fraction" yunits="fraction"/>
      <screenXY x="0" y="1" xunits="fraction" yunits="fraction"/>
    </ScreenOverlay>
  </Document>

But I cannot load this into Google Earth properly.  Frankly, I'm not sure if the KML file is wrong or if my workflow is wrong.  What I do know is that if I run the model scripts, R writes a file to the cwd and then looks for it elsewhere:

> plotKML(eberg["CLYMHT_A"])
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Writing to KML...
Closing  eberg__CLYMHT_A__.kml
The file /Users/flip/Google Drive/2017.summer/CMIP5/RCP_data/Bulk/R26_spatial/Drive/2017.summer/CMIP5/RCP_data/Bulk/R26_spatial/eberg__CLYMHT_A__.kml does not exist.

Moreover, on my Mac, this bring sup filehelper in the Google folder (I created because it seems to need one) and doesn't know that the KML file should be ahndled by Google Earth (possibly because that'd a chrome app now).


texts:
map = new OpenLayers.Map('map', options);
    // BASELAYER

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

    var google_satellite = new OpenLayers.Layer.Google(
        "Google Satellite", {
            type: google.maps.MapTypeId.SATELLITE,
            numZoomLevels: 20

        }
    );
    // OVERLAYS
    tiled = new OpenLayers.Layer.WMS(
    "Kaski",
    wmsUrl,
    {
        LAYERS: 'KaskiGrup',
        transparent: true,
        visibility: true,

    },
    {
        transitionEffect: 'none',
        isBaseLayer: false,

    }
);

    // NOKTASAL LAYERS
    wfs_layer = new OpenLayers.Layer.Vector("Point",
        {
            /*style: {
                externalGraphic: 'http://www.openlayers.org/dev/img/marker.png',
                graphicWidth: 21,
                graphicHeight: 25,
                graphicYOffset: -24
            },*/
            styleMap: new OpenLayers.StyleMap({
                "default": new OpenLayers.Style({
                    pointRadius: 10, // sized according to type attribute
                    fillColor: "#00FFFF",
                    strokeColor: "#33CCCC",
                    strokeWidth: 3,
                    graphicZIndex: 1
                }),
            }),

            strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],

            protocol: new OpenLayers.Protocol.WFS({
                version: version,
                // loading data through localhost url path
                url: wfsUrl,
                featureNS: featureNS,
                visibility: true,

                // layer name
                featureType: featureType,
                // geometry column name
                internalProjection: WGS84,
                externalProjection: WGS84_google_mercator,
                geometryName: geometryName,
                schema: schema,

            })
        });

I want to add google map on the background and overlay my wms and wfs layers but my wfs layers is not at correct location, also wms layers is at correct location. As you can see, i add google prejection on map option where is my mistake ? Thanks


texts: Computing large distance matrix with pgrouting?

texts:
The last weeks I have tried to make a Distance Matrix with pgrouting.

On small numbers I get a result. For instance < 100 will be 1 sec. < 1000 aprox 15 min. and 2000 aprox 20 min.  If I am trying to go for larger  numbers, for instance 4000, my computer crashes.... I have tried to change shared buffers, temp buffers, but that seems to make the system more unstable.

I am using postgreSQL 9.6, pgrouting 2.4.1   on a WIN 7 machine with 8GB RAM.
On Pgadmin query tool.

SELECT start_vid, end_vid, agg_cost FROM pgr_dijkstraCostMatrix(
    'SELECT gid AS id, source, target, trunc(length_m) AS cost, trunc(length_m) AS reverse_cost FROM public.ways',
    (SELECT array_agg(id) FROM  public.ways_vertices_pgr WHERE id < 100)
);

Does any one have some advice?

I would like a matrix to try before I would go for a compromise like Euclidean.


texts:
I'm installing QGIS 1.7.4 in windows 7 & having problems with QSpatialite. I get 3 messages...

Message #1 This DB a valid SQlite DB, but not a valid SpatiaLite DB. QSpatiaLite will now try to Create a valid SpatiaLite DB

then Message #2 Error. Try to execute by yourself: SELECT InitSpatialMetadata()

then Message #3 You are using SpatiaLite V2.3 ( standart with QGIS<1.7). Please note that QspatiaLite full potential require SpatiaLite V2.4. Please update your SpatiaLite version, or move to QGIS 1.7

What I've tried

Executing the SELECT InitSpatialMetadata() command from within the dialog
I'm running QGIS 1.7.4
The QSpatialite dialog confirms the error message because in the bottom left corner it says I'm running Spatialite V2.3.1. Thinking my problem relates Spatialite version, I ran the OSGeo4W setup again. However, the setup says I'm running spatialite 3.0.1-2 and sqlite 3.7.10-2.

The only thing I can think of is that spatialite v3 is not being picked up because the QSpatialite plugin installed v2.3.1 of spatialite somewhere. Is that possible? Any ideas?

Thanks

Ian


texts:
Is there a option / way in ModelBuilder to move the files I am processing in a Folder to a specific other folder (like the input folder - and done folder)?

I now Iterate files in a Folder, make a table of it,  and append to a geodatabase in a server.

So I don't need to move them everytime myself when it's all done.


texts: How can I draw a polygon with a fixed area in ArcMap

texts: Is there a way to draw a polygon with a fixed area in ArcMap?  For example, I'm trying to create a 40,000 foot boundary around a feature, but the boundary does not have to be any specific length or width.


texts: Accessing KML for Google Map API using Google Fusion Table

texts:
I am trying to host my KML file on fusion table so that I have dynamic control over my KML file and I am trying to access same as in public URI as an parameter in the Google Map Java Script .

I am getting lost on how to do it need some help on that .


texts: What can be used to drape an High res Image over a DEM and get high quality results?

texts:
I have a High resolution Georeferenced Orthomosaic,  as well as a DEM of that area. The pixel size of both is 20 cm. The data is in Tiff format.

I tried to drape the Image on the DEM in ArcScene, but the output doesn't have the resolution of the input Image.

You can see that the image appears blocky and pixelated. This pit must be about 40 meters in diameter, and looks much better in the orthomosaic.

Are there any other software which can be used to drape a georeferenced image over a DEM ?

While I'll prefer a free software, I don't mind a proprietary source, paid software.


texts: Flatten broad area through road fill connecting upstream and downstream

texts:
What are some recommended approaches for flattening a broad area through a road, as a means of connecting the upstream and downstream parts of a stream bed?

I am not trying to connect these areas for flow accumulation or stream network computation, but to model the area as a dam within HEC-RAS - so single cell connections, that would work for watershed / flow accumulation focused analyses, will not work here.

The area I wish to flatten, preferably smoothed from upstream edge to downstream edge, is shown below in blue. Aerial imagery with LiDAR overlaid.


texts: Renaming pre-existing table Field name using PyQGIS

texts:
I am using PyQGIS with QGIS 2.84.

This is the code I have used:

layer = iface.activeLayer()

parking_layer_field_names = [str(i.name()) for i in point_layer3.pendingFields()]

layer.startEditing()
for key, field in enumerate(layer.pendingFields()):
    if str(field.name()) != parking_layer_field_names [key]:
        layer.renameAttribute(key, 'new')
iface.vectorLayerTools().stopEditing(layer)

The error thrown by this code:

'QgsVectorLayer' object has no attribute 'renameAttribute' See log for more details


texts: Postgis : Simplify a path using ST_Simplify

texts:
I have a table with latitude, longitude representing paths.
I want to simplify the paths using Postgis ST_Simplify

Here is a sample set of 89 GPS points i try to simplify : https://gist.github.com/vdaubry/ab57027d8100495bdf261478808a8a50

I expect the result to be a smaller list of GPS points representing a similar shaped path.

When i run ST_Simplify on this set i only get 2 points as a result.

Can you help me debug this query ?

SELECT ST_X((result.points).geom) AS latitude, ST_Y((result.points).geom) AS longitude
FROM (
  SELECT ST_DumpPoints(simplified_path.simplifiedLine) AS points
  FROM (
    SELECT ST_Simplify(ST_MakeLine(path.point), 3) as simplifiedLine
        FROM (
            SELECT ST_SetSRID(ST_MakePoint(latitude, longitude),4326) AS point
            FROM tour_points
            where tour_id=369
            ORDER BY passing_time ASC
       ) as path
    ) simplified_path
) result

P.S: I tried playing with the tolerance parameter but i always end up with only 2 points.


texts:
I'm filtering a geojson layer by user inputs,and works fine for the first query, and works fine for the first query, the problem is with the second so on. starting from the second query brings me the last layer plus the new one I already tried map.removeLayer(layer) any ideas on how can I get only the new filter?,

function query(data) {
    var input=null;
    var user = null;
    var input = document.getElementById("userInput").value;

    departamentos.then(function(data) {
    var alldeptos = L.geoJson(data);

    user = L.geoJson(data, {

        filter: function (feature, layer) {
            //console.log(feature)
            return feature.properties.nom_dpto == input;
        },

    });
    });
    $("#allbus").click(function() {
        map.addLayer(alldeptos)
        map.removeLayer(user)
    });
    $("#form").click(function() {
        map.removeLayer(alldeptos)
        map.removeLayer(user)
        map.addLayer(user)
    });
   });
}


texts: Programmatic use of Geonetwork 3.1 metadata insert service

texts:
I'm trying to use the xml.metadata.insert service on an installation of Geonetwork 3.1. I've had success putting the parameters (a snippet in the case of data) directly in the URL in a web-browser, such as the correct response is returned (uuid and id of the created resource) and a metadata record is created.

However, trying to emulate the same thing in either python (through the "requests" package) or in postman, results in a strange behaviour of a 200 status code, but html being returned (below) instead and the resource not being created.

<!DOCTYPE html
  SYSTEM "html">

        <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
        <title>ICES Metadata Catalogue - International Council for the Exploration of the Sea (ICES)</title>
        <meta charset="utf-8">
        <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
        <meta name="apple-mobile-web-app-capable" content="yes">
        <meta name="description" content="">
        <meta name="keywords" content="">
        <link rel="icon" sizes="16x16 32x32 48x48" type="image/png" href="../../images/logos/favicon.png">
        <link href="rss.search?sortBy=changeDate" rel="alternate" type="application/rss+xml" title="ICES Metadata Catalogue - International Council for the Exploration of the Sea (ICES)">
        <link href="portal.opensearch" rel="search" type="application/opensearchdescription+xml" title="ICES Metadata Catalogue - International Council for the Exploration of the Sea (ICES)">
        <link href="/geonetwork/static/gn_login_default.css" rel="stylesheet" media="screen">
        <link href="/geonetwork/static/bootstrap-table.min.css" rel="stylesheet" media="screen">
        <link href="/geonetwork/static/ng-skos.css" rel="stylesheet" media="screen">
        <link href="/geonetwork/static/srv_custom_style.css" rel="stylesheet" media="screen">

Is there anything that I'm doing incorrectly to receive this response?


texts: QGIS convert long integers stored as string to integer

texts:
I have long integers (11 characters) stored as text in a shapefile that I need to convert to integers, but when I attempt to convert them in the QGIS field calculator, I get the following error:

Eval Error:
Cannot convert '26001980100' to int

The value in the error code above is the value from the first row in the field. I though this might be due to a field width constraint in shapefiles for integer columns, so I then imported the shapefile into a spatialite database and tried again without success. In both cases I was trying to create a new field populated by the converted values. How can I convert this field to integer?

The problem appears to goes deeper: after creating a new field and populating it with the integer values via the spatialite command line tool, QGIS does not then "see" the new field in the updated table. Am I out of luck with long integers in QGIS?

QGIS 2.12.3

Edit: I'm not sure whether this is a duplicate because the error appears specific to cases of converting strings having a length greater than 10 to integer. I tried the answer supplied here to convert to integer without luck. However, I was able to use that answer to convert the field to REAL, which is adequate for my purposes. The field in question is a foreign key, and the tables containing the primary key all store the values as INT, so I wanted to be consistent, but I can still join INT to REAL.


texts: Error Loading spatialite extension in .NET Project using System.Data.SQLite

texts:
I am trying to use Spatialite in a .NET application using the
System.Data.SQLite provider.

I have tried to load 'libspatialite-4.dll' but an error happens when I
call function "conn.EnableExtensions(true);conn.LoadExtension("libspatialite-4.dll");"

The specific error is "Attempted to read or write protected memory.
This is often an indication that other memory is corrupt."

I am using sqlite-netFx40-setup-x64-2010-1.0.92.0,libspatialite-4.1.1, VS2010 SP1, .net Framework 4.0, Windows7 64bit,


texts:
I am the developper of gisgraphy. I plan to update my reverse geocoding service to house numbers (it was only for streets).

Understanding what reverse geocoding is seems very simple and wikipedia says:

  Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name

But; imagine I am on a highway (let's say street A) and there is a house in 5 kilometers (house A), an hotel ( hotel A) in 1.5 kilometers on that road, and there is an street intersection (street B) in 2 kilometers where there is a house (house B) in 1 kilometers (so house B is at 3km from where i am).

If I reverse geocode my current position, what do I should  expect:

The street A (cause house A is too far and reverse geocoding is only about address, not POI ),
The hotel A, cause it is the nearest place on the road where I am, and geocoding is not only about address but also about POI)
The house A cause it is the nearest house on the road I am
The house B (cause it is the nearest house, even if it is not on my road)

What are the usages?

I am sure it depends on what you want, but what reverse geocoder usually do ?

I already have a street search API by position and a find nearby (POI), should I simply search the nearest 'feature' (street, house, poi, whatever...) or limit geocoding to address only ?


texts:
I am trying to update the symbology of a layer with another layer with a python script. I have written a code for this.
But I am getting the following error message :

ERROR 000732: Input Layer: Dataset C:\Temp\Assessment\work\code Output does not exist or is not supported
ERROR 000968: The symbol layer does not match the input layer

Here is a part my code that is experiencing this.

# Apply symbology from a layer

SymbologyLayer = arcpy.mapping.Layer("C:\Temp\Assessment\mod05AssessedEx\Shadow      Line Symbols.lyr")
Layers = arcpy.mapping.ListLayers(SymbologyLayer)

inputLayer = arcpy.MakeFeatureLayer_management(RasterIntoPolygon,"newLayer")
saveHere = r"C:\Temp\Assessment\work\code Output"
arcpy.SaveToLayerFile_management(inputLayer,saveHere)
print arcpy.Describe("C:\Temp\Assessment\mod05AssessedEx\Shadow Line Symbols.lyr")
DEMwithSymbol = arcpy.ApplySymbologyFromLayer_management(saveHere,Layers[2])
arcpy.RefreshActiveView()
arcpy.RefreshTOC
print "Applied symbology"

Can anyone please give any suggestion???


texts: QGIS, making a topo map from an excel file

texts:
I will appreciate your help. I am new to QGIS, this question might be common but trough me search I did not find the answer I needed hence my question here.

How do I make a topo map ( raster) from the excel file (X,Y UTM), attached is an a image since I did not have a better idea how to upload some of the data.
Thanks


texts: Referencing ArcMap layer to get spatial reference

texts:
I'm currently using ArcGIS 10.2.1. In previous custom geoprocessing tools, I've used the following code to create a feature class in-memory, and set the spatial reference using an existing layer in ArcMap:

polygon = arcpy.GetParameterAsText(0)

mem_point = arcpy.CreateFeatureclass_management("in_memory", "mem_point", "POINT", "", "DISABLED", "DISABLED", polygon)

Where polygon is a layer in the ArcMap TOC. This has always worked in the past, but I just went to use one of the tools, and I got the following error:

Traceback (most recent call last):   File "", line 23, in <module>
     mem_point = arcpy.CreateFeatureclass_management("in_memory", mem_name, "POINT", "", "DISABLED", "DISABLED", polygon)   File "",
line 1800, in CreateFeatureclass
     raise e ExecuteError: ERROR 000622: Failed to execute (Create Feature Class). Parameters are not valid. ERROR 000628: Cannot set
input into parameter spatial_reference.

Did this just recently change?


texts: QGIS 2.14.3 constantly crashing when digitizing

texts:
I'm using QGIS 2.14.3 Essen on Linux Mint 17.2 Rafaela. I'm digitizing either into .OGR vector file or a Spatialite database. Same happens with both.

I can be digitizing a string, I complete it with right mouse button, I then add in the attributes (automatically created unique id) and save. I then save the element. So, typically I'll then select another string nearby - with 'Select Feature by Area or Single Click' and then I crash straight out. If I haven't saved the line string I've lost it of course.

The message I get from the terminal is:

  "ERROR 4: Permission denied qgis.bin: geos_ts_c.cpp:3657: int
  GEOSCoordSeq_setOrdinate_r(GEOSContextHandle_t,
  geos::geom::CoordinateSequence*, unsigned int, unsigned int, double):
  Assertion `0 != cs' failed."


texts: Filtering feature from Atlas with PyQGIS

texts:
I need to print a layout that shows densities from data in a layer with a Graduated Symbol that has to be reclassified with quantiles based on the data of a field from the layer; this field is shown by atlas.  For example, county=1 and density=1 to 10, and county=2 density= 5 to 20, and so on.
My code exports several times the files and overwrites them (file1 county=1 densities 1 to 10; file2 county 2 densities 2 to 5; file1 county 2 densities 1 to 10; file2 county 2 densities 2 to 5,...).
My code up to now:
projectInstance= QgsProject.instance()
projectLayoutManager = projectInstance.layoutManager()
for comp in projectLayoutManager.printLayouts()
    if comp.atlas().enabled():
        comp.atlas().beginRender()
        nbobj = comp.atlas().updateFeatures()
        comp.atlas().first()
        for i in range(0, comp.atlas().count ()):
            print(str(comp.atlas().currentFeatureNumber()))
            active_layer.setSubsetString('"COUNTY"={}'.format(i+1))
            active_layer.renderer().updateClasses(active_layer,QgsGraduatedSymbolRenderer.Quantile,10)
            for legendLyr in iface.mapCanvas().layers():
                renderer = legendLyr.renderer()
                myRenderer = renderer.clone()
                legendLyr.setRenderer(myRenderer)
                legendLyr.triggerRepaint()
            comp.atlas().setFilterExpression('"COUNTY"={}'.format(i))
            comp.atlas().filterExpression()
            exporter =  QgsLayoutExporter(comp.atlas().layout () )
            exporter.exportToImage(comp.atlas(), "d:/path/", 'png',QgsLayoutExporter.ImageExportSettings())
        comp.atlas().endRender()

How do I export only the desired atlas feature in every iteration ("i" feature, so file1 county=1 density= 1 to 10 and file2 county=2 density=2 to 5,...), instead of all the features for every iteration?


texts:
I have seen on several maps in the case where a polygon overlaps (or maybe underlaps) the data frame, which is distributed to the margins, and text is turned on, that the text can print outside the data frame in the margin area.

It shows this on both the pdf export and the print straight to a plotter.
Has anyone seen this same behavior?


texts: Error adding feature to centroids in QGIS

texts: I want to map forest stands. I have the compartments as centroids and I am trying to add a feature by snapping to them and creating polygons. Each time I click, however, it displays a box titled "Compart...", with an text box after the letters DN. It doesn't matter what I enter, I can't seem to create a polygon.


texts: Export a KEA segmented file as shp or raster, using segment ID numbers

texts:
I have a KEA dataset of which the RAT has numerous columns. There are 11919 segments.  I want to export this dataset either as shp or raster keeping the segment ID numbers, i.e. I want a shapefile or raster with only the segment numbers connected to the segments.

How do I do this?

I tried rastergis.exportGDALCols.Image in Python using RSGisLib but since the segment ID number column doesn't have a column name, that didn't work.
I was able to export the RAT as a csv in Python with an additional column that contains segment numbers, but I need a shp or raster dataset to give me the spatial information of the segments.


texts: Style according to column name in CartoCSS for a combined set of columns

texts:
I am trying to create a label for a points layer in Carto using names of three different column names.
I figured that the following code works to get the correct behaviour.
#layer::labels  {
  text-name: [name] + " - " + [address]  ;

  text-face-name: 'DejaVu Sans Book';
  text-size: 20;
  text-fill: #FFFFFF;
  text-label-position-tolerance: 0;
  text-halo-radius: 1;
  text-halo-fill: #385871;
  text-dy: -10;
  text-allow-overlap: true;
  text-placement: point;
  text-placement-type: dummy;
}

However, I also would like to give the columns a different color. For example 'name' would have lets say red and 'address' will have orange.
I tried the following.
  #layer::labels  {
  text-name: [name] + " - " + [address]  ;
  text-face-name: 'DejaVu Sans Book';
  text-size: 20;
  text-fill: #FFFFFF;
  text-label-position-tolerance: 0;
  text-halo-radius: 1;
  [name]{text-halo-fill: #385871;}
  [address]{text-halo-fill: #fcde9c;}
  text-dy: -10;
  text-allow-overlap: true;
  text-placement: point;
  text-placement-type: dummy;
}

But the above gives syntax errors. How can I do this?