texts: Applying graduated symbol renderer only on some features using PyQGIS

texts:
I have set up a graduate symbol renderer (QGIS cookbook - Graduated Symbol Renderer), and I am now trying to make it conditionally applied to some features only, depending on one (or more) attributes. For example, only trying to apply the render to features having X value and not display the other ones.

Is QgsRuleBasedRenderer to be used or should I go with creating a virtual layer with only the required features following my requirements?

from qgis.PyQt import QtGui

file = 'C:/Users/Administrator/Desktop/20200123_T00_QGIS.txt'
uri = "file:///" + file + "?delimiter=|%5Ct%7C&skipLines=27&skipEmptyFields=Yes&detectTypes=yes&xField=Longitude&yField=Latitude&crs=EPSG:4326&spatialIndex=no&subsetIndex=no&watchFile=no"

Hsep = QgsVectorLayer(uri, 'Hsep', "delimitedtext")
QgsProject.instance().addMapLayer(Hsep)
layer = Hsep
targetField = 'H-sep'

rangeList =[]
opacity = 1

# define value ranges
minVal = 0.1
maxVal = 0.3
# range label
lab1 = '0.1 - 0.3'
# color (yellow)
rangeColor = QtGui.QColor('light green')
# create symbol and set properties
symbol1 = QgsSymbol.defaultSymbol(layer.geometryType())
symbol1.setColor(rangeColor)
symbol1.setOpacity(opacity)
symbol1.symbolLayer(0).setStrokeStyle(Qt.PenStyle(Qt.NoPen))
#create range and append to rangeList
range1 = QgsRendererRange(minVal, maxVal, symbol1, lab1)
rangeList.append(range1)

# define value ranges
minVal = 0.3
maxVal = 0.5
# range label
lab2 = '0.3 / 0.5'
# color (yellow)
rangeColor = QtGui.QColor('red')
# create symbol and set properties
symbol2 = QgsSymbol.defaultSymbol(layer.geometryType())
symbol2.setColor(rangeColor)
symbol2.setOpacity(opacity)
symbol2.symbolLayer(0).setStrokeStyle(Qt.PenStyle(Qt.NoPen))
#create range and append to rangeList
range2 = QgsRendererRange(minVal, maxVal, symbol2, lab2)
rangeList.append(range2)

# create the renderer
groupRenderer = QgsGraduatedSymbolRenderer('', rangeList)
groupRenderer.setMode(QgsGraduatedSymbolRenderer.EqualInterval)
groupRenderer.setClassAttribute(targetField)
# apply renderer to layer
layer.setRenderer(groupRenderer)

# add to QGIS interface
QgsProject.instance().addMapLayer(layer)


texts: Basemap tiles are disjointed due to bad projection

texts:
I am working within a corporate environment that is using ArcGIS Server (v10.2.2). I am new to trying to work with ArcGIS API for JavaScript (v3.9) to create custom map applications. I can't seem to get a basic map up and running. For example, when I write this code:

<!DOCTYPE html>

<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no">
<title></title>
<link rel="stylesheet" href="MYDOMAINHERE.COM/3.11/esri/css/esri.css">

I get the basemap map tiles all disjointed like so (re-created in Photoshop for illustration purposes):

Is it possible that the coordinates are wrong?

Our server is using ArcGIS with 4326(WGS84).


texts: Querying for all rows that have the pattern of a constant value in one field but different values in another field using ArcMap

texts:
In ArcGIS (ArcMap 10.8 or Pro 2.6) I have a table of 4000 buildings (rows). Each row/building belongs to one of fifty groups and thus has a group ID of 1-50 in the field GROUP_ID. Each row also has a GROUP_STATUS field which can have the value of "fully open" or "partly open".
Some groups have conflicting statuses. For example, there are multiple rows in Group_ID = 10, 14, or 40 which have values of both "partly open" and "fully open" in GROUP_STATUS. For these, and all other groups that have the same issue, I need to change the status to only "partly open". I need to find out for which GROUP_IDs this is the case so I can update the status accordingly.
Without doing a frequency table - only a selection or definition query - how do I identify all GROUP_IDs that have different values for GROUP_STATUS in the table?
I would like to find rows where GROUP_ID = a single value and GROUP_STATUS = more than one value across the rows, along the lines of "select by GROUP_ID where GROUP_ID is consistent but GROUP_ID is not consistent."


texts: Batch GDAL_Calc multiple band calculation

texts:
I have a question concerning gdal_calc.

I´m running windows 7 64b and am working from the OSGeo4W shell.
I have multispectral images that need to be prepared for analysis but they are not saved in a single band per wavelength format. Instead each wavelength is saved in a tif 8b format and this needs to be converted to a 16b tif format with a single band. The data is stored in the green and blue bands and a formula is available to calculate the new band. Images are acquired with Tetracam Mini-MCA system and more details can be found here.

So my question is: how to use gdal_calc to batch process images?

The only manual solution I found is separating the bands using gdal_translate and repeating the calculations for each image separately.


texts: Knowing Previous Command Item C# in ArcObjects?

texts:
I am creating a tool for ArcGIS Desktop that will execute a command when I click a place on the map.

The Code looks like this:

public class TileLasLoader_Tool : ESRI.ArcGIS.Desktop.AddIns.Tool
{
    private static EsriTools.Forms.TileLasLoader_Form _form;
    ICommandItem _previousCommand = null;

    public TileLasLoader_Tool()
    {
        _form = EsriTools.Forms.TileLasLoader_Form.instance;
    }

    protected override void OnUpdate()
    {
    }

    protected override void OnMouseDown(MouseEventArgs arg)
    {
        base.OnMouseDown(arg);
        if (arg.Button == MouseButtons.Left)
        {
            _form.LoadLasFile(arg.X, arg.Y);
            ArcMap.Application.CurrentTool = _previousCommand;
        }
    }

My Only problem is getting the tool to revert back to the previous "Command Item" that it was before I used the tool. If it was not using a tool I would just make it null so it reverts back to the windows mouse cursor.


texts: only one popup displayed on the openlayers

texts: The issue is that i only one popup on the map is dispalyed. I can see all the points but popup can see only on the last point of the list. Please check my code http://jsfiddle.net/3HSjY/24/ I'm stack with rhis problem from Sunday and can not find any info that could solve my problem on the internet :(


texts: Select Multiple feature using select interaction

texts: How to select multiple feature using select interaction in openlayers 3?Is it possible. Multiple feature can be selected using drag-box. But I have needed using select interaction.


texts: Can OpenLayers show map in projection other than wgs84 or UTM, for example Albers projection?

texts: I've managed to show map in both utm or wgs84 in OpenLayers but Albers projection seems not possible as it does not have bound? Not sure about this, might be what I'm doing is wrong.


texts:
I'm new to GIS, having only dabbled with ArcGIS for Desktop for about 4 weeks in 2007. The business I am at wants to target writing complete solutions that involve GIS as it would benefit their hardware sales. They do not know which market they want to target but want me to pick the GIS software to use, yet I am very green to GIS markets. All i know if target platforms are windows, windows mobile, android and web applications. The idea of mobile applications means data connection is not 100% so internet based services may need to be ruled out, i.e. GoogleMaps.

I was struggling to work out how to phrase this as a question, given the huge range of GIS products available, Esri's ArcGIS for Desktop, CadCorp, MapInfo, GoogleMap/Earth, qGIS, Manifold, openJump, openstreetmap etc.

To avoid the what do I buy question, I instead ask the question why do you use Esri ArcGIS. From my experience its slow to navigate around maps and quite verbose to program in causing quite a lead time in development. Esri is also the most expensive product on the market, and i believe also the most popular. But why do people continue to use it, when other systems like CadCorp have the same functionality if not more, seem to be more responsive, can import Esri data and cost less.

I feel that I must be overlooking something here. I know client demands in some industries like defense, want Esri but why is it still so popular, given I often hear how other products are betters (Is this a Windows vs. Mac example on market share vs. quality)

P.s. This question has been spawned from this one on StackOverflow.


texts: PostGIS check if multiline string is included in polygon

texts:
I want to return all records that the multiline-string is included in the polygon. In the photo below, it should return two records.
Green line (streets) is contained inside two neighbourhoods (white line)
I use ST_Contains(a."NeighbourhoodWatchGeometry", b.geom); in the WHERE clause. It doesn't return any record for that street.
FROM "NeighbourhoodWatch" a,
    gis_streets b
  WHERE
st_contains(a."NeighbourhoodWatchGeometry", b.geom);


texts: How do I calculate a weighted centroid with respect to a raster, e.g. population weighted centroid

texts: In ArcGIS 10.0, how do I calculate the weighted centroid of a collection of polygons, where I want the weight to be a raster layer. In particular I am using a raster of population counts and would like to generate the population weighted centroid as a result.


texts:
I tiled my mxd document using the Multithreaded MapCruncher sample located here: http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//000100000mm8000000

It appears to work, the output directory has different level directories, L0, L1... which contain .bundle files.

I don't know how to connect to these .bundle files.  I have looked at the DynamicCacheLayerManager, but that needs to be initialized with a layer file and the .bundle files are created with an MXD file which contains 1 or more layers (in my case 2 layers).  I tried it but it didn't work.   How can I connect to .bundle files?  I want my application to read these files and load the images instead of re-tiling the loaded maps.

I'm using ArcEngine 10 C# VS2010 Thanks in advance.

So I guess this page:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//000100000mmm000000
says the following,

Once the cache is finished processing, open Engine or Desktop, navigate to the cache folder, and add the cache as a raster.

I tried using the "Open" button and the "Add Data" button and navigating to my .bundle files, but neither of these recognizes .bundles files.


texts: Raster interpolation with gdalwarp, and replicating with python API

texts:
My problem is similar to this: Replicating result of gdalwarp using gdal Python bindings. I want to resample a raster using bilinear interpolation, which I've done successfully with gdalwarp by adding height and width specifications for the output file:

gdalwarp -ts w_output h_output -r bilinear -co 'TFW=YES' -of VRT input.tif output.vrt

I'm trying to replicate the above result using the Python API. Can it be done using gdal.AutoCreateWarpedVRT?

src_ds = gdal.Open(input.tif)
src_wkt = None
dst_wkt = ????
resampling = gdal.GRA_Bilinear
output_vrt = gdal.AutoCreateWarpedVRT(src_ds, src_wkt, dst_wkt, resampling)

I'm guessing I have to do something clever with the dst_wkt in order to replicate the above gdalwarp result. I just want to increase the number of pixels by a factor of 10 compared to the input raster.

P.S. If it's better to just use gdalwarp, that's good to know. I have no problem doing that.


texts: Not getting correct result from for loop

texts:
The for loop below works without error, but the result is get in the table of my database is wrong. Its as if the for loop skips right over the first 3 if and elif statements and goes to the else. The table I am using has values 1 through 5 in each row of the WWFLOWDATA column and the data type is a number. Does anyone have any idea why its not giving me the output I am looking for? Maybe an issue with the access database?

import os
import arcpy
from arcpy import env

arcpy.env.workspace = r"C:\Users\sutton\Documents\Python\PythonExcercise\DuncanvilleData.mdb"

try:

        #create new field to hold new values
        arcpy.AddField_management("Wastewater_Lines","PipeAssesment","TEXT","50")
        print("Rehab Field added to Wastewater_Lines")
        with arcpy.da.UpdateCursor("Wastewater_Lines",["WWFLOWDATA","PipeAssesment"]) as cursor:
            cntr = 1
            for row in cursor:
                #update the indicator field
                if row[0] == '5':
                    row[1] = "Critical"
                elif row[0] == '4':
                    row[1] = "Monitor"
                elif row[0] <= '3':
                    row[1] = "Unknown"
                else:
                    row[1] = "Unknown"
                cursor.updateRow(row)

except Exception as e:
    print(e.message)

StatusofScript = "Script executed"
print StatusofScript


texts: How to add Marker on map click using OpenLayers 6?

texts:
I want to implement adding icon/marker on map click. It should work like this (Leaflet example). https://codesandbox.io/s/add-onclick-on-mapcontainer-in-react-leaflet-v-3x-4r7tc
I am using OpenLayers 6.5 and React (java script tip would be fine enough thought). My current code looks like this. My idea was to use Draw Api to style the point to Icon manner but it failed to work. Line and Polygons are styled correctly thought. Is there any efficient way to do this or any workarounds?
Additional it would be ideal to swap current default point style to my Icon (following cursor movement).
const drawSource = new VectorSource({ wrapX: false });

class PublicMap extends Component {
  constructor(props) {
    super(props);

    this.state = {
      center: [0, 0],
      zoom: 1,
      draw: null,
      drawVectorLayer: new VectorLayer({
        source: drawSource,
        style: new Style({
          fill: new Fill({
            color: `rgba(4, 149, 106, 1)`,
          }),
          stroke: new Stroke({ color: "blue" }),
          image: new Icon({
            anchor: [16, 48],
            imgSize: [32, 48],
            src: "http://openlayers.org/en/v3.20.1/examples/data/icon.png",
            zIndex: "10",
          }),
        }),
      }),
    };

    this.olmap = new Map({
      layers: [
        new Tile({
          source: new OSM(),
        }),
        this.state.drawVectorLayer,
      ],
      view: new View({
        center: [0, 0],
        zoom: 2,
      }),
      target: null,
    });
  }

  componentDidMount() {
    this.olmap.setTarget("map");
  }

  addInteraction(e) {
    this.olmap.removeInteraction(this.state.draw);

    if (e.target.value === "None") {
      this.setState({ draw: null });
      return;
    }
    this.setState(
      { draw: new Draw({ source: drawSource, type: e.target.value }) },
      () => {
        this.olmap.addInteraction(this.state.draw);
      }
    );
  }

  undoLastPoint() {
    this.state.draw.removeLastPoint();
  }

  render() {
    return (
      <>

        <form class="form-inline">
          Geometry type: &nbsp;
          <select id="type" onChange={this.addInteraction.bind(this)}>
            <option value="None">None</option>
            <option value="Point">Point</option>
            <option value="Polygon">Polygon</option>
            <option value="Circle">Circle</option>
          </select>
          <input
            type="button"
            value="Undo last drawn point while drawing polygon. "
            id="undo"
            onClick={this.undoLastPoint.bind(this)}
          />
        </form>
      </>
    );
  }
}

export default PublicMap;


texts: Measuring minimum distance between objects in two layers using ArcPy?

texts:
I am working on an alien species analysis looking at the spread of a species from the mainland to islands in an archipelago. I am doing the analysis using arcpy and I have succeeded in isolating risk zone islands based on how far the alien species can swim. Yup, the critters can swim some 2.000 meters.

With my islands isolated using repetitive buffer analysis I now need to be able to measure the minimum distance between one of the islands in the risk-zone and the mainland. So what I am looking for is something like this:

arcpy.distancebetweenlayers(risk_islands,mainland)

The closest I can get is the point distance tool. Since it only works with points I am out of luck. "Exploding" the features into points could be an option, but I would rather not have to go that way.

Is there some functionality I have overlooked which can support me on this task?
Using other libraries than arcpy might be option. Which libraries could give me a quick fix on this one?

For the record I am aware of a similar question with response (Get shortest distance between two geometries in ArcPy) about a year ago, but do not think it sufficiently answered the question.


texts: How to highlight routes using a plugin like ATLAS/EASY PRINT

texts:
First I'm new to GIS and this forum.
Secondly sorry my bad English.

I use Atlas to auto generate my maps.  Eevery "route" has a coverage.
At times the "routes" are too close, so the lines cross and it is impossible to make out a specific "route". I really need some advice on this please.

Atlas already jumps from "route" to "route" take the .jpeg images for me.

I tried to use EASY PRINT but was unsuccessful.

The black is what I get by default. The red is what I would like to show for that "route" with the highlighted "route" name (red added with paint).

So, in short words, when that map gets generated for that "route" it should be highlighted and the rest of the "routes" must stay black.  This needs to be done automatically: I have too many maps to do this manually.

I'm sorry if this topic has been covered but I did search and can't find help (probably due to lack of correct search, and poor English).  I hope you guys understand and can help me as soon as possible. Thank you !


texts: Displaying information in infowindow using KMLLayer and kmz file?

texts:
I added the kmz layer to my map

www.spc.noaa.gov/products/md/ActiveMD.kmz

However when I click on it, all I get in the infowindow is "MD 1558". (This may change for you since it is updated in near realtime)

If I download the kmz file and rename to a kml and open in Notepad I see the following code below.

Is it possible to include the http://www.spc.noaa.gov/products/md/MD1558.kmz link in the infowindow? If not, is it possible to configure an infoWindow with the KMLLayer so that I could add my own link?

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document>
<Folder>
<name>SPC Active MDs as of Fri Aug 15 18:11:02 UTC 2014</name>
<open>1</open>
<description>
<![CDATA[
Visit http://www.spc.noaa.gov/products/md/ for more info.
]]>
</description>
<NetworkLink>
<open>1</open>
<name>MD 1558 Active Till 1915 UTC</name>
<visibility>1</visibility>
<refreshVisibility>1</refreshVisibility>
<description>
<![CDATA[
MD 1558 Active Till 1915 UTC.<br />
Visit http://www.spc.noaa.gov/products/md/md1558.html for more info.
]]>
</description>
<Link>
<href>http://www.spc.noaa.gov/products/md/MD1558.kmz</href>
<refreshMode>onInterval</refreshMode>
<refreshInterval>60</refreshInterval>
</Link>
</NetworkLink>
</Folder>
</Document>
</kml>


texts:
I have the exact same problem as explain in this post. I'm trying to find the shortest path between each building on a map (shp multipolygon) to any transformation plant (shp multipoints) by following
a network (shp multipolyline).

As suggested in the answer given to the OP, I was thinking of using pgr_dijkstra() but I think it won't work. The thing is, if I'm not mistaken, pgr_dijkstra() return the shortest path between 2 points. In this case, we do not know which electric transformation plan we're looking to reach, the goal is to find the shortest path to any transformation plant on the network.

In order to use pgr_dijkstra(), we would have to first assign an electric transformation plant to each building (minimum radial distance around it for example) but there could be situations in which the closest transformation plant in radial distance is not at all the closest by the network.

Does anyone knows if there is a way, with pgRouting or something else, to find the shortest path to any of the electrical transformation plant on the network (I've already managed to connect each building to the closest line on the network in order to be able to look for the shortest path between points which already are on the network, cf. this post).


texts: Updating and interacting with map layers through a windows form

texts:
Using arcobjects I am attempting to connect a windows form with a checkedlistbox, to the layers that I have on my map in ArcGIS. When I first create my form, the list of checkboxes will automatically populate itself with the name of all the layers, and the checkbox next to the layer will be checked if the layer is set to be visible in the arcmap, and unchecked if it is not.

I want to make my form so that if the user checks/unchecks a box in the form, the corresponding layer will check/uncheck on my arcmap (and vice versa). I am having trouble getting the value (checked or unchecked) of a layer from the checkedlistbox once I try to change it on the form, and I also don't know how to update the windows form once I do so. If I use form.show(), it comes up with an entirely new map, and I would have to repopulate the checkbox and do other things. I wanted a way to simply update the values on the existing form that I had. Someone suggested using form.update(), but that isn't doing what I want it to.


texts: How to convert elevation data from Google Elevation API to NAVD 88?

texts:
I want to convert the elevation data obtained from Google elevation API to NAVD 88. I understand that Google uses WGS84 with EGM96 Geoid.
Do you know what I should do to convert the elevation data? Also, is there any python library that could help me with the conversion?


texts:
I am trying to map a WFS layer with data originally coming from a PostGIS table, and using a OpenLayers 3. At some point in this process, I believe the longitude and latitude coordinates are getting switched. Here is the relevant bit of code, but I don't think this is the problem.
 var planningAppsSource = new ol.source.GeoJSON({
      'projection': map.getView().getProjection(),
      'url': 'http://localhost:8080/geoserver/wfs?DATASTORE=Framework&service=wfs&version=1.3.0&request=GetFeature&typeName=POSTGIS:equatortable&outputFormat=application/JSON&maxFeatures=5'
 });

var planningAppsLayer = new ol.layer.Vector({
   source: planningAppsSource,
   style: new ol.style.Style({
      image: new ol.style.Icon(({
        anchor: [5, 5],
        anchorXUnits: 'fraction',
        anchorYUnits: 'pixels',
        src: 'marker-icon.png'
      }))
   })
});

map.addLayer(planningAppsLayer);

None of the points showed up, so I went and entered the output from the WFS call into geojson.io, and points that were supposed to be in Indiana were in Antarctica. Switching the x and y in geojson.io got the points to where they should have been. To double check my hypotheses I created a PostGIS table with the point (0,0), (5,5) and (-5,5). The first two were where you would expect them to be, while the last was not. The odd things is that the WFS is the only place where I have this problem. I can call the same exact data using WMS, and it maps fine using openlayers or on the geoserver previews. I can get the WFS and view it using QGIS and the map looks right. Any idea where the problem is?
My best guess is that the problem arises because I have my data is in EPSG:4326 projection which takes lat, lon while the protocol for WFS 1.3 is lon, lat. Am I on the right track for understanding what is wrong?
I can probably fix the problem by creating 'false' points, and having
everything plotted back properly, but I was wondering if there was a less absurd solution to this problem.


texts: Add to display Feature Class (GDB) with variable name

texts:
I have a model in modelbuilder and I'm trying to display the output with the same name as the input (using Parse Path Name) which something like "C12345_BIRDS" so I want to add "_FINAL" to this. If I set a concrete name everything is displayed properly including the symbology but if I set the output name as %Value%_FINAL the feature class is created as well but not displayed. The parameters and tick-boxes are correctly set.

I have read somewhere this is a ESRI bug and it has to be solved by using Code Block box and arcpy.mapping in Calculate Field, that I used right after I created the output I'm interested in.

I dont know what to write in there to force my model to display a feature class with variable name.

Any idea?


texts: How does Virtual File System (VSI) File handling work in ArcGIS Pro?

texts:
There are a few references to GDAL VSI urls in the ESRI documentation, for example here but no documentation on how it's implemented or used.
For instance, adding  /vsizip//vsicurl/https://download.geofabrik.de/europe/andorra-latest-free.shp.zip into the add Data From Path does not work in ArcGIS Pro, but I can add this string into the Add Vector Layer in QGIS and it loads without any issue


texts:
Does anyone know how to perform a overpass-api request like this?

"Show me all banks in Germany"?

I could just show all 'amenity = bank' in a specific bounding-box. But Germany (or most of all countries) doesn't fit very well in such a bbox.

Is there any way to do this? Have you an example? Server performance or runtime is secondary since we could host a own server or something like that :)


texts: How to simplify a multiline geometry in QGIS?

texts: I'm trying to simplify my shapefile, which uses a lot of vertices to represent curvatures, by using the Generalizer plugin. But no matter what threshold I use for Douglas algorithm, it still includes almost all of those closely placed vertices, which I can see are less than 1 feet away from each other. Any hints of what I'm doing wrong?


texts: Selective attributes not displaying in GeoCortex Essential Manager

texts:
I'm having trouble with Geocortex Essential Manager displaying the full set of attributes for a particular feature. Specifically, I have a parcels dataset and 4 related tables. I recently added a new subdivision into the base and filled out the related tables for the parcels.

Normally, I just need to reconcile my edits against the default and not touch GeoCortex Essential Manager and the changes will appear on the webmap. However this time only the geometry edits (parcel polygons and address points) went through and appeared in the webmap. The data in the tables was not viewable in the webmap.

I went back and checked the realtes in Arc and everything worked great, I've made sure all the changes were reconciled back to the default version on the database and I've been thru Geocortex and reviewed all the options and nothing seems to be off. Every other parcel in the base still has the full range of attributes associated with it.

Is there a lag time between edits being reconciled in Arc and the changes appearing in the webmap?


texts:
I have the following Model which is working fine (1st image).

I'm running a Network Analyst Location-Allocation. I want to iterate both Facilities and Demand Locations.
So I'm calling a sub model (2nd image) to iterate the demand locations, which is stored as differents shapefiles within one folder.

However, when the model finishes it only used one shapefile from that folder and not all 4 shapefiles that it contains.

Can you figure out what I'm missing?


texts: Import QGIS Python modules in Python command line or standalone script

texts:
I am running a Windows 7 machine with Python version 2.7.8 (32 bit).  I installed QGIS  2.4 (Chugiak) from OSGeo4W and simply want to be able to import the QGIS Python modules (e.g. qgis.core) into the previously installed Python, in much the same way one can do so with arcpy.

I have Googled and set various folders in my PATH and PYTHONPATH environment variables but nothing seems to be able to let me say:

import qgis.core

I just get the ol' module not found error.  Has anyone else had success?

Thanks.


texts: Which EPSG setting when georeferencing German maps from 1880s?

texts:
I'm trying to georeference a map from 1881 that was created by the Preußische Landesaufnahme. When using Gauß-Krüger (EPSG:31467), my map section always ends up way too far in the east. I've found out by now that they used another meridian back in the days, the Ferro Meridian.

Is there a preset EPSG in QGIS for the Preußische Landesaufnahme to compensate the differences to todays coordinate systems?


texts: Problems Integrating Sub-Model into Master Model

texts:
I am attempting to integrate a iterative model into a master model to merge all of the outputs of the iterative model into one feature class. However, when I attempt to set the output of the sub-model (iterative) as the input for the merge tool, it will not allow it. I have attached two pictures of each model. Any suggestions?

I am using ArcMap 10.

Thanks,

Eric


texts:
When I execute the Erase Tool in ArcGIS it does not erase all of the features in the input class. The Input and Erase features are both polygons. I have attempted multiple times, hoping it was just a fluke.

Here is what I have tried so far:

1) I have selected out the specific features from the input that aren't being erased and attempted to run the Erase with those, but the output shapefile ends up being completely empty.

2) Using the Clip tool instead. The effects are the same.

3) Made sure when I run the Erase that no random features are accidentally selected.

EDIT:

4) Repair Geometry of input.

5) Rebuild Spatial Index of input.

6) Making sure input featues (erase featues) are not multipart via a)Dissolve and b) Multipart to Singlepart Tool. The latter produces a shapefile that is missing the 3 polygons that are causing trouble during the Erase function.

Does anyone have experience about how to fix this problem?

Input features:

Output feature. The bottom half of the polygon is where the trouble is:


texts: Parameters for python script: how to provide a user interface that lists a tool's available value options?

texts:
I have several scripts that I have working from inside PythonWin. I am now experimenting with running these from inside ArcMap 10.0 so I can have a user define the parameters of the script. An important part of designing this popup for the user is choosing a "data type" for each parameter. However, I find ESRI's information about data types for parameters of little use for figuring out which to choose to get different results for user interaction.

For instance, say I have a script uses the dissolve tool in some way and I would like to have the user choose what type of "statistics_fields" they want in the final resulting shp file that the script returns (there are nine choises: FIRST, LAST...etc). Specifically, I would like to know how to choose the correct data type for this parameter so that the final user interface:

a) Provides a list of all nine choices available and allows the user to pick up to 9 variables that they would like to be included in the attribute table of the resulting shp file.

b) Provides a list of all nine choices as a pull down menu and allows the user to only pick one that ends up in the attribute table of the resulting shp file.*

How can I accomplish a) and b)? Perhaps I need to do more than just choosing the right data type?


texts: UpdateLayer Function Resulting in AssertionError

texts:
I am providing an alternative to the accepted answer posted here. I pose as new question as the script I use is resulting in error and I've been unable to resolve.

Esri help documentation states that the UpdateLayer function, when altering all of a feature classes symbology properties, may be used as an alternative to the RemoveLayer and AddLayer functions. The script below is intended to meet the same end as the question posted above, but it uses this UpdateLayer function, which seems like it would be a more elegant solution than using the RemoveLayer and AddLayer functions.

When I run this script:

import os
from os import listdir
from os.path import isfile, join
import arcpy

# Creates list of MXDs
mxdList = [f for f in listdir(r'C:\MXDs')
             if isfile(join(r'C:\MXDs', f))]

for doc in mxdList:
    mxd = arcpy.mapping.MapDocument(os.path.join(r'C:\MXDs', doc))

    # FCs for which there are .lyr files
    fcList = ['FC1', 'FC2', '...']

    # Iterate through dataframes
    for df in arcpy.mapping.ListDataFrames(mxd)[0]:
        print df
        workSpaceType = "FILEGDB_WORKSPACE"

        for item in fcList:
            for lyr in arcpy.mapping.ListLayers(mxd, item):
                desc = arcpy.Describe(lyr)
                geometryType = desc.shapeType

                # Use geometry to match lyr file name
                if geometryType == 'Point':
                    # Access and define the lyr file
                    newLyr = arcpy.mapping.Layer(os.path.join(r'C:\lyrFiles', item + '_P.lyr'))
                    print 'Layer Object:   ' + str(newLyr)
                    arcpy.mapping.UpdateLayer(df, lyr, newLyr, False)
                    print 'Updated:  ' + str(newLyr)
                    del newLyr

                elif geometryType == 'Polyline':
                    # Access and define the lyr file
                    newLyr = arcpy.mapping.Layer(os.path.join(r'C:\lyrFiles', item + '_L.lyr'))
                    print 'Layer Object:   ' + str(newLyr)
                    arcpy.mapping.UpdateLayer(df, lyr, newLyr, False)
                    print 'Updated:  ' + str(newLyr)
                    del newLyr

                elif geometryType == 'Polygon':
                    # Access and define the lyr file
                    newLyr = arcpy.mapping.Layer(os.path.join(r'C:\lyrFiles', item + '_A.lyr'))
                    print 'Layer Object:   ' + str(newLyr)
                    arcpy.mapping.UpdateLayer(df, lyr, newLyr, False)
                    print 'Updated:  ' + str(newLyr)
                    del newLyr

    mxd.saveACopy(r'C:\MXDsUpdated' + doc)
    del mxd

it results in error:

assert isinstance(data_frame, DataFrame)
    AssertionError

for line 58, for the UpdateLayer function for the Polygon geometry, the geometry of the first item in the fcList list for the first MXD that is accessed.

I have ensured I have defined the Layer as described here. Also, inserting the print df statement within the ListDataFrames for loop prints feature class names and not dataframe names. Any ideas?


texts: Attribute Lost with CopyFeatures_management()

texts:
I have a Network Analysis Layer and I would like to export one of its classes as a .lyr file. Doing this in the GUI is easy enough, just right click and export as a layer.

However, in arcpy and ModelBuilder, I have been trying to use CopyFeature replicate this process. It works for the most part, but the exported layer loses the ObjectID attribute (the temporary file before FeaturetoPoint also is missing ObjectID), which is retained when I do this through the GUI.

Here is the code exported from ModelBuilder:

dir = "C:\\Users\\Ariel\\Dropbox\\Maps\\Geo_Data\\"
Origins = "OD Cost Matrix\\Origins"
v170614_OD_Cost_Matrix_Final_lyr__2_ = dir +
"170614_OD_Cost_Matrix_Final.lyr"
temp_origin_shp = dir + "temp_origin.shp"
v170623_Origins_shp = dir + "170623_Origins.shp"

# Process: Copy Features
tempEnvironment0 = arcpy.env.extent
arcpy.env.extent = 170614_OD_Cost_Matrix_Final.lyr (2)
arcpy.CopyFeatures_management(Origins, temp_origin_shp, "", "0", "0", "0")
arcpy.env.extent = tempEnvironment0

# Process: Feature To Point
arcpy.FeatureToPoint_management(temp_origin_shp, v170623_Origins_shp,
"CENTROID")

I could probably try to export the attribute table and then join it back in in python, but might anyone know why CopyFeatures fails in this manner?


texts:
I have a set of polygon centerlines created using voronoi polygons. The method seems to work well for what I want to accomplish, but these 'tail' remnants are throwing off successive analysis (attached an image). Is there a way to remove them from the lines automatically? Note: I would prefer if this could be done with QGIS tools, but can create a python script if necessary.

I would be open to switching methods for obtaining polygon centerlines. I know there is a 'Polygon Centerline' tool in ArcMap, but I don't have access to it.


texts: GeoServer gives No XSDElementDeclaration found Exception

texts:
I have added the following WFS source to the Stores in GeoServer. http://geoservices.rijkswaterstaat.nl/verkeersscheidingsstelsel_noordzee?service=wfs&version=1.1.0&request=GetCapabilities

This went well, and I can now see the list of supported feature types in the web interface, when trying to add a new layer. However, when selecting one of these feature types, I get an error. The feature type I have tried to add in this case is named nzvss_ank (named wfs_nzvss_ank in GeoServer). The relevant part of the stack trace is:

Caused by: java.lang.RuntimeException: org.geotools.data.DataSourceException: No XSDElementDeclaration found for {http://www.opengis.net/wfs}nzvss_ank
    at org.geotools.data.store.ContentFeatureSource.getAbsoluteSchema(ContentFeatureSource.java:342)
    at org.geotools.data.store.ContentFeatureSource.getSchema(ContentFeatureSource.java:309)
    at org.geotools.data.store.ContentFeatureSource.getSchema(ContentFeatureSource.java:109)
    at org.geoserver.catalog.CatalogBuilder.buildFeatureType(CatalogBuilder.java:348)
    at org.geoserver.catalog.CatalogBuilder.buildFeatureType(CatalogBuilder.java:330)
    at org.geoserver.web.data.layer.NewLayerPage.buildLayerInfo(NewLayerPage.java:323)
    ... 97 more

The full stack trace can be found in this Gist.

The same exception occurs when trying to add other feature types, and when using other WFS links from the same provider.

Other information: GeoServer version is 2.6.0, charset set for the XML documents is UTF-8, feature buffer size is set to 10.

What can I do to fix this issue?


texts: Making QGIS layer update from changed data source

texts:
I'm trying to get layers to update automatically when their data source changes. I'm using R to write a shapefile with an attribute, and colouring according to that attribute in QGIS.

I want to write a new shapefile with different attribute values, and have the Qgis map colours update. Step 1 is triggering that process, step 2 is making the layer reload from the modified shapefile. Its step 2 I am worrying about here.

Other questions/mailing list chatter mentions using triggerRepaint on the layer - that doesn't work. Other suggestions include setCacheImage(None) and again that doesn't work. The layer does update eventually, but I really can't see the logic, and sometimes it happens by surprise after I've done nothing. Or maybe I did something two minutes ago.

The one reproducible way of getting it to update is to duplicate the layer from the legend menu - the duplicate always gets its data from the current shapefile, and the original layer updates itself too! So there must be some way of doing it.

I think it was working better in 2.8, but this is 2.10 so maybe there's a new bug somewhere.

Related, but doesn't work for me in 2.10:

How to automatically reload raster layers if source is changed in QGIS?

Other things I've tried:

layer.dataProvider().dataChanged.emit() - worked once, then not again on the same layer

I think I've tracked down why duplicating the layer works - if I create a new throwaway layer based on the updated layer and then call .triggerRepaint() on the updated layer, it updates on the map canvas:

QgsVectorLayer( layer.source(), "layer copy", layer.providerType() )
layer.triggerRepaint()

If I use a different layer source it doesn't work, so it seems to be if you create a layer object based on the same layer source...

A quick test just now with a raster layer (from a GeoTIFF), and just calling rlayer.triggerRepaint() seems to reliably update the view of the raster in the map canvas.


texts: Convert polyline shapefile to points 1m apart and add xy coords in R

texts:
I have a table with xy coordinates in R and the following packages loaded: spatstat, raster, rgdal, sp and zoo. I do know how to use this to make a point shapefile, a line shapefile that goes through these points and a smoothed line shapefile using rollmean on my original xy coordinates. My workflow it is something similar to this:

Table with xy coordinates
Build a spatial object from table
save the shapefile shape points giving the name of the spatial object, the directory where the file will be saved and the name of the file without extension.
Use same table with xy coordinates to do a polyline shapefile from a SpatialLinesDataFrame
Using rollmean do a smoothing of the xy coordinates, and save it as a smoothed polyline shapefile.

Ok, now I have all the objects I want, including SpatialLines and SpatialLinesDataFrame for the smoothed coordinates.

Now I want to do the following 2 main transformations using R (I know how to do that if I use ArcGIS but I don't want to go back and forth from R to ArcGIS):

a.  Transform the smoothed shapefile (or any of the above "smoothed" objects) into points at 1 m apart, and add to the shapefile an attribute table with xy coordinates for each point along the smoothed line

b.  Transform the smoothed shapefile in p points equally distributed along the smoothed line (that means for me equal distance – although I don't know what that distance is) and add to the attribute table xy coordinates for each point.

# make shapefiles - point, line, smoothed line
#tabl is a table with the first 2 columns x and y coordinates which is what is used
#to actually do the shapefile
#rast - the raster from which we take the projection and the resolution
# smooth.win = how big you want the smoothing window, the real smooth param is actually the
# no. of cells you want to smooth over
#dir.name = path of the directory where i want to save the shape files.
#Need to use double back slash when writing the directory path name.
# file name – a name that will be added to the name of the raster that is analysed

make.shape <- function(tabl, dir.name, file.name, smooth.win, rast){

require(spatstat)
require(raster)
require(rgdal)
require(sp)
require(zoo)

crs <- projection(rast)
res.x<-round(xres(rast),1)
res.y <- round(yres(rast),1)
rast.name <- names(rast)
file.name <- paste(rast.name, "_", file.name, sep = "")

if(is.data.frame(tabl)) tabl <- tabl else tabl <- as.data.frame(tabl)

# build spatial objects that can be saved as ESRI shapefiles
shape.sp <- SpatialPointsDataFrame(tabl[,1:2], tabl, proj4string=CRS(crs))

# save the shapefile shape points giving the name of the spatial object,
#the directory where the file will be saved and the name of the file
#without extension.

writeOGR(shape.sp, dir.name, paste(file.name, "_points", sep = ""),
   driver = "ESRI Shapefile")

# PolyLine shapefile

xy <- cbind(tabl[,1], tabl[,2])
cl <- Line(xy)
cl1 <- Lines(list(cl), ID = 1)
row.names(tabl) = seq(1, dim(tabl)[1])

shape.l <- SpatialLines(list(cl1), proj4string = CRS(crs))
shape.ln <- SpatialLinesDataFrame(shape.l, tabl)

writeOGR(shape.ln, dir.name, paste(file.name, "_line", sep = ""),
     driver = "ESRI Shapefile")

# do the smoothing of the shape line
rmy <- rollmean(tabl[,2], smooth.win)
rmx <- rollmean(tabl[,1], smooth.win)

x1y1 <- matrix(c(rmx,rmy), ncol=2, nrow = length(rmy), byrow = FALSE)
cl1a <- Line(x1y1)
cl1b <- Lines(list(cl1a), ID = 1)
row.names(x1y1) <- seq(1, dim(x1y1)[1])
shape.l1 <- SpatialLines(list(cl1b), proj4string = CRS(crs))
x1y1 <- as.data.frame(x1y1)
shape.ln1 <- SpatialLinesDataFrame(shape.l1, x1y1)

writeOGR(shape.ln1, dir.name, paste(file.name, "_smooth_line_",
     round(smooth.win*res.x,0), "m", sep = ""), driver = "ESRI Shapefile")
}


texts: Correct CRS for viewing LandScan data in QGIS?

texts:
I think my problem is simple, but I'm still stumped. I'm using QGIS 2.14.3 (having recently switched from Arc). I have downloaded LandScan data for a small area (northern Virginia) and am attempting to bring some XY data onto my map.

LandScan is in WGS1984 and, as instructed in the landscan documentation, I set my project CRS to be WGS 84. But for the life of me, I cannot get my XY data (nor my TIGER lines) to show up in the right spot relative to my landscan raster. LandScan documentation is very clear that the LandScan data should not be reprojected, lest you will lose data, and I know my TIGER files are properly projected, so I'm quite confused as to how to proceed. I think I might be able to figure this out in Arc, but as I don't currently have access to Arc, I have to do it in QGIS.

I'm grateful for any help!


texts:
There seems to be a newly formed issue in displaying (?) XYZ tiles in QField.

In preparing a project in QGIS (desktop) I have for perhaps the last year successfully viewed XYZ tiles in QField after having added tiles in to a QGIS project via the directions outlined in the answer of @MatthiasKuhn (~3rd answer from the top) to this post: Adding Basemaps from Google or Bing in QGIS?

(See also: https://www.lutraconsulting.co.uk/blog/2016/10/26/qgis-xyz-tile-wmts-preview/)

However, recently the XYZ tile layers, while working fine on the desktop, do not either display or copy over in or into QField. (they do seem to be listed in the QField application's legend, though they are "greyed out")

I have usually used the QConsolidate plugin to transfer my project, but have now tested with the QFieldSync plugin also.

If the methods described in the linked answer above are no longer a successful way to access XYZ tiles (such as google satellite) in QField, what is a current way to access XYZ tiles in a QField project?

[could the web addresses for the tiles just have changed?]

Using QField 10.4 - Kesch on Android 5.1.1 from a project prepared on QGIS 2.18.15 (Win 7).


texts: ArcMap/ArcObjects - Drawing graphics onscreen: Explanation or tutorial for Draw/Refresh/PartialRefresh methods?

texts:
I find myself woefully lacking in understanding of when, why, and how (which type) to do screen draw/refresh/partial refresh.  For instance, in a standalone Python script outside the application boundary, I have a function that is trying to draw graphic elements in the data view.  Specifically, the elements will be a TextElement on top of a RectangleElement.  The function loads the two elements into a GroupElement, which is in turn added to the ActiveView's GraphicsContainer for drawing on the screen.  (As far as I can tell the graphic elements are valid.)  I am really not understanding what Draw/Refresh/ParitalRefresh methods to put where.  If anyone can give/point to a good explanation or tutorial about screen drawing/refreshing I would appreciate it.  Below is a bit of selected code from my function.  I have a PartialRefresh in the last line, but I don't know if that is correct or not.

pGroupElement = NewObj(esriCarto.GroupElement, esriCarto.IGroupElement3)
pTextElement = NewObj(esriCarto.TextElement, esriCarto.ITextElement)
pRectElement = NewObj(esriCarto.RectangleElement, esriCarto.IElement)
.
.
(configure TextElement)  [tedious code removed]
.
.
(configure RectElement)
.
.
ptxtElem = CType(pTextElement, esriCarto.IElement)

pGroupElement.AddElement(ptxtElem)
pGroupElement.AddElement(pRectElement)

pContainer = mxDoc.ActiveView.GraphicsContainer

pGrpElem = CType(pGroupElement, esriCarto.IElement)

pContainer.AddElement(pGrpElem, 0)

mxDoc.ActiveView.PartialRefresh(esriCarto.esriViewGraphics, None, None)

Added:  So, is creating objects with the following code - as I do now - not creating them in process?

def NewObj(MyClass, MyInterface):
    """Creates a new comtypes POINTER object where\n\
    MyClass is the class to be instantiated,\n\
    MyInterface is the interface to be assigned"""
    from comtypes.client import CreateObject
    try:
        ptr = CreateObject(MyClass, interface=MyInterface)
        return ptr
    except:
        return None

If not, I gather I must use IObjectFactory.  Assuming that my current objects are not in process, I must say I haven't seen any great performance hit, but I guess I'm just asking for trouble forcing lots of interprocess calls.

Thank you, all, for the help.  I shall experiment when next I get the chance and let you know the results.


texts:
Is there a way to save the zoom, tilt, and orientation settings of a 3D map view in QGIS?

QGIS version 3.10.

I have a 3D map view that I need to save as Image file (e.g. using the "Save as Image..." function).

I need to do this many times (once for each layer) and over different time periods.

I need the view geometry (zoom, tilt, and orientation) to be always exactly the same each time I save the image. Eventually, it would be nice to have more geometries saved and to be able to quickly switch from one to another.

Is there a way to save the zoom, tilt, and orientation settings so to reuse them?


texts: Request data from OpenLayers to a standalone GeoWebCache

texts:
I'm actually using OpenLayers and GeoServer to make a GIS application. For the moment, GeoWebCache is directly integrated into GeoServer, but for my future needs, I'll need to have it as a stand alone.
I've dowloaded it and deployed it, but now I can't make it work.

But I basically want to do is, for the moment, use it as a "firewall".
Typically, in OpenLayers I have such requests :

"http://129.182.247.82:8080/geoserver/wms,"

to retrieve WMS data.

What I want to do, is to change this URL to make it point to the GeoWebCache stand-alone, and then the GeoWebCache stand-alone will call the right layers into GeoServer.

I followed this documentation : http://geowebcache.org/docs/current/quickstart/index.html

But I can't make it work. I had no geowebcache.xml to test layer so I created one from github and put it into WEB-INF/classes, but it still does not work.
Inside my geowebcache-core-context.xml I've added lines like this :

<bean id="gwcWMSConfig" class="org.geowebcache.config.GetCapabilitiesConfiguration">
...
<constructor-arg value="http://http://129.182.247.82:8080/geoserver/ows?service=WMS&amp;request=GetCapabilities&amp;version=1.1.0">

and this

<bean id="gwcTLDispatcher" class="org.geowebcache.layer.TileLayerDispatcher">
....
<constructor-arg ref="gwcGridSetBroker"/>
<constructor-arg ref="gwcWMSConfig"/>
<constructor-arg ref="gwcXmlConfig"/>
</bean>

But still nothing is working, and I'm quite lost :'( !
Any help would be gladly appreciated :-)!


texts: Reference raster data outside of GeoServer data directory

texts:
is there a possibility to point GeoServer to raster data stored neither in a database nor in the data directory of GeoServer?

My problem is that IT keeps me from direct access to the GeoServer machine and I only have access via a terminal machine. That keeps me from storing the data in the data directory.


texts: Unmanaged Raster Catalog Relative Paths?

texts:
Is it possible to store relative paths in an unmanaged raster catalog?

What I want to do:

C:\Folder\File.gdb
    UnmanagedRasterCatalog1
        ..\..\Data\MyRaster\Raster.File

C:\Data
    MyRaster
        Raster.File

I see this feature request which leads me to believe it's not possible:

Relative paths in unmanaged raster catalogs - http://ideas.arcgis.com/ideaView?id=08730000000brlG

When I use the Repair Raster Catalog Paths tool it accepts a relative path as the new path parameter but when I export the paths using the Export Raster Catalog Paths tool the paths are the full path. So it looks like it is resolving the full path and putting that in instead of the relative path.

Repair Raster Catalog Paths tool

Export Raster Catalog Paths output


texts:
I have a problem converting a shapefile to raster.

I created a point-shapefile using a text-table with some xyz-coordinates. Now I want to convert it into a tiff, where all the areas without values (everything except the points) get the value 0.

When I use the tool "rasterize" I can only decide for one attribute (x, y or z) with a black box without any values as output. What am I doing wrong?


texts:
I am working on modify function for an Openlayers project. I define a vector layer for modify objects. (line, polygon etc.) But then I take an error like below while adding vector layer to map object. Why am i getting this?

Error

Uncaught TypeError: layer.getLayerStatesArray is not a function
at ol-debug.js:17395
at ol.Collection.forEach (ol-debug.js:8289)
at ol.layer.Group.getLayerStatesArray (ol-debug.js:17394)
at ol.Map.ol.PluggableMap.renderFrame_ (ol-debug.js:18670)
at ol.Map.<anonymous> (ol-debug.js:17624)

Source Code (core)

        init: function (Map) {
            var self = this;
            self.Map = Map;
            self.VectorSource = new ol.source.Vector({ wrapX: false });
            self.Map.olMap.addLayer(self.VectorSource);
            self.initialized = true;
        },

Source Code (UI)

var geo = (new ol.format.GeoJSON()).readFeatures(e, { featureProjection: 'EPSG:3857' })[0];
Max.Map.Draw.VectorSource.addFeature(geo);

Referenced example


texts: Add "Position and additional attributes" in table

texts:
Is there any way to add to a shape attribute table "position and additional attributes" provided while tracking with GPS?
I mean, it could be useful to add automatically in the table, for each feature collected attributes like 'altitude', 'speed', 'direction', 'quality' and so on.


texts: Best Practices for Tracking Feature Changes and Supporting Documentation

texts:
In my organization, it is very important that we have documentation from an authoritative source to support every feature that we add to our data and every edit that we make to each feature. Currently, our organization keeps an excel file on which we track every change we make and include links and file paths to the documentation that supports the change. Lately, these excel spreadsheets have become a bit cumbersome.

I am curious if anyone else deals with this level of documentation, and how they do it. I know about metadata and the track edits tool, but neither seems to provide the level of detail that we need to maintain data lineage and documentation for every single feature.

Is there a tool or a preferred method for dynamically linking each feature in a feature class to detailed information about it's lineage and documentation?

Thanks!


texts: Generate points in multiple features taking minimum distance into account

texts:
I have a function which creates Wind Turbines represented as points. Essentially, it uses the code from the Random points inside polygons (fixed) tool albeit with some slight changes.

The goal is to create random points inside polygons taking into account the minimum distance specified. This works very well especially with polygons which are not close to another (e.g. a single polygon):

However, if the polygon is close to or adjoined to another polygon (e.g. as shown below), the points from each polygon can be within the minimum distance as shown in red:

How could I change the code so that those points in red are not close to another from a nearby polygon?

Ideally, I would like multiple points to be replaced by a single point:

Here is the code which can be reproduced in the Python Console, a polygon layer must be selected with a relevant CRS before running the function:

import random
from PyQt4.QtCore import QVariant

def checkMinDistance(point, index, distance, points):
    if distance == 0:
        return True
    neighbors = index.nearestNeighbor(point, 1)
    if len(neighbors) == 0:
        return True
    if neighbors[0] in points:
        np = points[neighbors[0]]
        if np.sqrDist(point) < (distance * distance):
            return False
    return True

def generate_wind_turbines(spacing):
    layer = iface.activeLayer()
    crs = layer.crs()
    # Memory layer
    memory_lyr = QgsVectorLayer("Point?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "Wind turbines for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    # Variables
    point_density = 0.0001
    fid = 1
    distance_area = QgsDistanceArea()
    # List of features
    fts = []
    # Create points
    for f in layer.getFeatures():
        fGeom = QgsGeometry(f.geometry())
        bbox = fGeom.boundingBox()
        pointCount = int(round(point_density * distance_area.measure(fGeom)))
        index = QgsSpatialIndex()
        points = dict()
        nPoints = 0
        fid += 1
        nIterations = 0
        maxIterations = pointCount * 200
        random.seed()
        while nIterations < maxIterations and nPoints < pointCount:
            rx = bbox.xMinimum() + bbox.width() * random.random()
            ry = bbox.yMinimum() + bbox.height() * random.random()
            pnt = QgsPoint(rx, ry)
            geom = QgsGeometry.fromPoint(pnt)
            if geom.within(fGeom) and checkMinDistance(pnt, index, spacing, points):
                f = QgsFeature(nPoints)
                f.setAttributes([fid])
                f.setGeometry(geom)
                fts.append(f)
                index.insertFeature(f)
                points[nPoints] = pnt
                nPoints += 1
            nIterations += 1
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

generate_wind_turbines(500)

Edit:

Dissolving and/or converting the polygons to singleparts does not seem to help much as the generated points still seem to fall within the minimum distance.

Tested on QGIS 2.18.3.


texts: Get relative north aligned point 2163 (US Equal Area)

texts:
What I'm trying to do is create a grid where the y axis is north aligned. I have a origin Lat/Lon point which is the bottom left corner of the grid and each cell should be 100 by 100 meters. So I projected my points to SRID 2163 (US Equal Area) and figured id use each the grid x,y values with the cell size relative to my origin point.

like this:

cell_bottom_left = { lon: origin.lon + (100 * x), lat: origin.lat + (100 * y)}
cell_top_left    = { lon: origin.lon + (100 * x), lat: origin.lat + (100 * y) + 100}
cell_bottom_right= { lon: origin.lon + (100 * x) + 100, lat: origin.lat + (100 * y)}
cell_top_right   = { lon: origin.lon + (100 * x) + 100, lat: origin.lat + (100 * y) + 100}

it worked fine, except the grid was tilted. I then found out that 2163 is not north aligned which explains what is happening here.

How would I go about creating a north aligned grid?


texts: windows function with points in a mixed (multipoint/point) table

texts:
I have a table comprising multipoints (tracks) and points (stops between tracks). there can be several consecutive multipoints.
now, i want to merge any points lying very closely to each other (spatially and temporally) and delete all multipoints lying between them.
if there would be only one multpoint and one point in turn, i could delete the correct multipoints with comparing the respective lead point with the respective lag point (window function), but in case of group_id's 434 - 437 there are three multipoints in a row. so, how to compare the next point with the previous point in each case, even if there are more than one multipoint between those points?


texts:
I'm a beginner in the world of geopackages.
I learnd from  here and here how to use gdal_translate (and gdaladdo) to import rasters in an existing geopackage or to create a new one. Works fine!
But how can I remove one raster that I don't want anymore in the dataset?

Edit after applying method described in user30184's answer: (I quote with copy-paste)
Step 1

$ ogrinfo -sql "drop table new_table" 20170831_Expertise_Elia_Gblx-Auvl.gpkg

INFO: Open of '20170831_Expertise_Elia_Gblx-Auvl.gpkg'
      using driver `GPKG' successful.

Step 2

$ ogrinfo -sql "delete from gpkg_contents where table_name='new_table'" 20170831_Expertise_Elia_Gblx-Auvl.gpkg

INFO: Open of '20170831_Expertise_Elia_Gblx-Auvl.gpkg'
      using driver `GPKG' successful.

Step 3

$ ogrinfo -sql "delete from gpkg_metadata_reference where table_name='new_table'" 20170831_Expertise_Elia_Gblx-Auvl.gpkg --config OGR_GPKG_FOREIGN_KEY_CHECK OFF

ERROR 1: bad result for PRAGMA foreign_key_check, got 11 rows, expected 0
ERROR 1: pragma foreign_key_check on '20170831_Expertise_Elia_Gblx-Auvl.gpkg' failed
ERROR 1: bad result for PRAGMA foreign_key_check, got 11 rows, expected 0
ERROR 1: pragma foreign_key_check on '20170831_Expertise_Elia_Gblx-Auvl.gpkg' failed
FAILURE:
Unable to open datasource `20170831_Expertise_Elia_Gblx-Auvl.gpkg' with the following drivers.
  -> PCIDSK
  -> netCDF
  -> etc...

=> At 3rd step it went wrong...


texts: Using ArcPy to Read from Excel Spreadsheet to Python Dictionary

texts:
How would I go about creating a script that reads a specified excel spreadsheet and places the values in a dictionary?

I have read about the xlrd library and how the library in such endeavors to read Excel spreadsheets in Python. However, I was hoping to strictly do this in the ArcPy module. Is this possible or would you be required to use the xlrd library?


texts: Preparing gridded data for GIS from Matlab

texts:
I'm looking for the best way to prepare raw data files for import into GIS. I am not a GIS user, but am analyzing data that will be passed to GIS users on our team.
I am performing an analysis in Matlab (but just as easily C++ or Python) for a grid of locations, resulting in a dataset with Latitude, Longitude, and Data (e.g. Hours of sunlight in a day). I could also describe the locations in X,Y coordinates for a stereographic projection with an arbitrary origin.
I'm dealing with polar data (e.g. origin 85S, 30E), so my Lat/Long is not square. X and Y coordinates are.
My intent is to write this data to files that can be directly imported into GIS (ArcMap). I believe this should be formatted as a raster layer. I'd like to include any required metadata to project it properly in GIS. I'm finding it quite difficult to find straightforward information on a format and syntax that would be appropriate.
Does anyone have any ideas or a good reference for data providers?


texts: Openlayers3 WMS GetFeatureInfo request only returns info for one feature

texts:
I have a point dataset in an Openlayers 3 map. The aim is to be able to click each point and get the information about that point displayed. The code below works just fine, it loads all the layers up and when I click a point, it retrieves some information. However it will only return information for one point, no matter where I click on the map. Can anyone see where I have gone wrong here?

var layers = [
  new ol.layer.Tile({
    source: new ol.source.MapQuest({layer: 'sat'})
  }),

  new ol.layer.Image({
    extent: [16935988.2390421, -3393028.96482803, 17102996.5023459, -3207079.63356002],
    source: new ol.source.ImageWMS({
      url: 'http://geonct.com:80/geoserver/NCTPROP/wms',
      params: {'LAYERS': 'NCTPROP:corr'},
      serverType: 'geoserver'
    })
  }),
  new ol.layer.Image({
    extent: [16935988.2390421, -3393028.96482803, 17102996.5023459, -3207079.63356002],
    source: new ol.source.ImageWMS({
      url: 'http://geonct.com:80/geoserver/NCTPROP/wms',
      params: {'LAYERS': 'NCTPROP:br_web'},
      serverType: 'geoserver'
    })
  }),
  new ol.layer.Tile({
      source: new ol.source.TileWMS({
          url: 'http://geonct.com:80/geoserver/NCTPROP/wms',
          params: {'LAYERS': 'NCTPROP:br_point'},
          serverType: 'geoserver'
        })
  })

];
function loadMap(){

var map = new ol.Map({
  layers: layers,
  target: 'map',
  view: new ol.View({
    center: [17031695, -3310901],
    zoom: 9
  })
});
var view = new ol.View({
  center: [0, 0],
  zoom: 1
});
map.on('singleclick', function(evt) {
  document.getElementById('info').innerHTML = '';
  var viewResolution = /** @type {number} */ (view.getResolution());
  var url = layers[3].getSource().getGetFeatureInfoUrl(
      evt.coordinate, viewResolution, 'EPSG:3857',
      {'INFO_FORMAT': 'text/html'});
  if (url) {
    document.getElementById('info').innerHTML =
        '<iframe seamless src="' + url + '"></iframe>';
  }
});
};
function layerSwap(evt){
   layers[evt.value].setVisible(evt.checked);
}


texts: Locking feature service without affecting its map service

texts: Is there a way to lock a feature service published by ArcGIS Server. For example: When I publish a service with feature access capability, I get two URLs; one for a map service, and the other for the feature service, but when going to the manager to lock the feature service, I can lock the map service, which locks both (the map and feature service). Is there any idea to lock the feature service only and keeping the map service unlocked?


texts: After joining shapefile with CSV file graduated symbols option unavailable in QGIS?

texts:
I have successfully joined data from a csv file to a shapefile in QGIS but for some reason I am having trouble symbolizing my data.

I want to use a graduated symbology however this option is not available for any of the columns for my joined table.

I have tried saving as a new layer but this does not change the issue.

Has anybody got any tips to get around this?


texts: I am using QGIS 3.2 and I need to create an area to form a polygon shapefile using latitude and longitude, not using Python as I am just a beginner. How to do this?


texts:
I'm modifying the Arrows plugin by Gregoire Piffault. It writes the output to a shapefile, but if the shapefile already exists it fails. He tries to remove the file using QgsMapLayerRegistry.instance().removeMapLayer then write a new one using QgsVectorFileWriter, which results in an empty file; if the file doesn't already exist it works fine.

I've tried simply deleting the file using os.remove after the call to removeMapLayer, but the file is still locked.

os.remove(self.conf["arrow_outshape"])
        WindowsError: [Error 32] The process cannot access the file because it is being used by another process: u'C:/GIS/SK/Saskatoon/mpArrows.shp

I think there are two options:

change the code substantially to allow for editing the existing file rather than creating a new one if it already is loaded
figure out how to completely unload the file after removing the map layer so that I can delete it


texts: Shapefile rotated after displayed on Bing Maps

texts:
I took a country shapefile from Thematic Mapping API | World Borders Dataset.

Then opened it in QGIS and saved needed layer as KML file and took coordinates for displaying polygons.

var country = new Microsoft.Maps.Polygon([[new MM.Location(35.360881805420263,46.286518096923942),...

In QGIS shapefile is showing correctly but on Bing Maps it is rotated and in another place.

What is the problem?


texts: Reading feature class attributes with or without quotes in ArcPy?

texts:
I'm working on a script that clips a batch of feature classes by a selection made in another feature class. The way it works is by reading the clip feature fields, then allowing the user to select the field they want to clip by. Then list of unique attributes fro the field is generated and the user can check multiple values to clip by. The multi value list is then passed to a SQL where clause that then gets incorporated into the SelectLayerByAttribute function.

The question that I've got doesn't really relate to the script (I've come up with some logic to get around the problem), but to how the attributes retrieved from the clip feature are formatted.

As I look at the values retrieved, some are formatted with single quotes around them ('value') and others are not (value). And both conditions can exist in the same feature class and field.

Can anyone tell me why this is the case and if there is something I should be doing with my data to eliminate this condition?

I've included the script that I've been using.  Note the if statement that includes .startswith and .endswith.  That's the logic I've been using to get around the problem.

import arcpy, os
workingGDB = arcpy.GetParameterAsText(0)
clipLyr = arcpy.GetParameterAsText(1)
clipField = arcpy.GetParameterAsText(2)
clipValue = arcpy.GetParameterAsText(3)
walk = arcpy.da.Walk(workingGDB, datatype="FeatureClass")
valueList = clipValue.split(";")
whereList = []
for v in valueList:
    if v.startswith("'") and v.endswith("'"):
        whereValue = '{0}'.format(v)
        whereList.append(whereValue)
    if not v.startswith("'") and not v.endswith("'"):
        whereValue = "'{0}'".format(v)
        whereList.append(whereValue)
arcpy.MakeFeatureLayer_management(clipLyr,"lyr")
whereClause = "{0} IN ({1})".format(clipField, ",".join(w for w in whereList))
arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION", whereClause)
records = arcpy.GetCount_management("lyr")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        filePath = os.path.join(dirpath, filename)
        outFile = filePath + "_clip"
        arcpy.Clip_analysis(filePath, "lyr", outFile)
        arcpy.Delete_management(filePath)
        arcpy.Rename_management(outFile, filePath)

Here is the validation code:

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""

    if self.params[1].value and self.params[2].value:
        self.params[3].filter.list = sorted({row[0] for row in arcpy.da.SearchCursor(self.params[1].value, self.params[2].value.value) if row[0]})

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return


texts: Get z value for a point of polygon using JTS library?

texts:
I have a polygon described with coordinates which have z value.

Is there any JTS utility that would allow me to get Z value for a coordinate which is inside the polygon?

I have a function for doing it.

Does it look reasonable?

private def getMinMaxEvelation(geometryWithZ:Geometry, someTwoDimGeometry: Geometry) = {
    val builder = new ConformingDelaunayTriangulationBuilder()
    builder.setSites(geometryWithZ)

    val zValues = someTwoDimGeometry.getCoordinates.toList.map{ coord =>
      builder.getSubdivision.locate(coord).dest().getZ
    }
    Elevation(minZ = zValues.min.toInt, maxZ = zValues.max.toInt)

  }


texts: How close project QGIS (* .qgs) without saving (PyQGIS)

texts:
I have several *.qgs projects and I want to open each project, add layer, export composer and close this project withot saving. Goto next project and repeat this operations.
But I did not find any close() method for QgsProject.
How I can close *.qgs without saving?
Something like this:

project = QgsProject.instance()
for prjFile in prjFiles:
   project.read(QFileInfo(prjFile))
   #do something with project
   #close project instance without saving and goto next project


texts: Eliminating empty values using arcade in ArcGIS Portal pop up?

texts:
I am attempting to create a conditional arcade expression that will list all of the fields that do not have a null value with an individual label and not show the field at all if the field is null. My current expression is below.

if (IsEmpty($feature.ABARTA_WI)){
return ""
} else {
return "Abarata: "+$feature.ABARTA_WI
}

This expression works in that it does not show the label or null value, but in the ArcGIS Portal pop up, it leaves an empty line for each one that is null, and I need it to not be there at all so the pop up is clean. As you can see below, this is creating empty lines in the popup.

I would prefer the popup to display something like below, eliminating the empty lines.

Chief: 0.28
Chief WI MKR: 0.28
Enerplus: 0.3


texts: Carto-CSS/ Mapnik XML Text background color

texts:
I want to add a background color to my labels. How do I do that with Carto-CSS in Tilemill or with Mapnik?

Carto-CSS

#marker_labels {
  text-name: "[name]";
  text-face-name:@sans;
  text-size: 20;
  text-fill: #FFFFFF;
  text-allow-overlap: true,
}

Mapnik XML

<TextSymbolizer face-name="Air Soft Cond Heavy" size="20" fill="#ffffff" allow-overlap="true" ><![CDATA[[name]]]></TextSymbolizer>

Here is a example on how the label should look like (ignore the line to the label):


texts: Using TiledMapServiceLayer and On Demand Tiling

texts:
I have a tiling services which generates tiles on demand but TiledMapServiceLayer does not use a callback for getTileUrl so how would one implement this?  I've tried using the refresh! method to no avail.  I've also tried using the undocumented img argument but didn't get far.  Thoughts?  Solutions?

Here's my code showing some of what I've tried:

getTileUrl: function ( level, row, col, esriImg )
{
    var that = this, imageNumber, folder, url;
    imageNumber = ( Math.pow( 2, level ) * row ) + col;
    folder = Math.floor( imageNumber / this.tileInfo.width );
    url = this.url.substitute( {
        x: col,
        y: row,
        z: level,
        folder: folder
    } );
    if ( h$.yes( mapVars.OnDemandTiles ) )
    {
        var img = dojo.create( "img", { src: url } );
        dojo.connect( img, "onerror", function ()
        {
            infor.gis.hansen.setMapEvent( "tilegeneration", {
                callback: function ()
                {
                    $.delay( 5000, that.id, function () {
                        that.refresh();
                    } );
                },
                custdata: url,
                extent: mapVars.FullExtent
            } );
        } );
        return url;
    }
    else
    {
        return url; // getRefreshURL(url);
    }
}


texts: Using fine-resolution samples to extract from coarser-resolution efficiently in R

texts:
I have a raster stack that I want to extract from:

raster resolution: 0.0417, 0.0417  (x, y).

My raster is coarser than my coordinates for extraction. I have finer coordinates whose resolution is ~ .000217, .000217. (My coordinates are coarser because I have other rasters with finer resolution). My coordinates aren't quite a rectangular grid (only am using coordinates in the form of a US state) but they increase incrementally in both LAT/LONG direction.

I've been using fast_extract from the library prioritizr & velox after converting my coordinates into SpatialPoints but I have too many coordinates and it's taking too long to go through each coordinate and extract from the raster.

I feel like the most effective way would be to extract at the resolution of the raster, and use that value for all MY coordinates that fall within that raster's pixel.

for example, let's take one pixel, which is from -74.0417 to -74.0834 (LONG) and 40 to 40.0417 (LAT). The value of this pixel would be for all my coordinates that fall in this range, ie. 40.00027 x 74.0417 / 40.00054 x 74.0417 / and so on.

I think I might want to use something like the aggregate function but I don't want to lose the number of my LatLon coordinates by making it coarser---I ultimately want to keep the same number of LatLon coordinates but reduce the number of extractions.

Here is some sample data

LatLon coordinates to extract:

temp_LAT<-seq(44.86,45.10,.0027) #89
temp_LONG<-seq(-80,-71.5,.0027) #3149
temp_LL<-data.frame(matrix(0,nrow = 280261, ncol =2))
colnames(temp_LL)<-c("LAT","LONG")
temp_LL[,"LAT"]<-rep(temp_LAT,each=length(temp_LONG))
temp_LL[,"LONG"]<-temp_LONG
pts<-SpatialPoints(temp_LL[,c("LONG","LAT")])
crs_prj<-"+proj=utm +zone=48 +datum=WGS84"
proj4string(pts)<-CRS(crs_prj)

raster to extract from:

library(raster)
x<-raster()
x <- raster(ncol=35, nrow=150, xmn= - 81.5, xmx= -71, ymn= 40, ymx=45)
res(x)<-.03333
res(x)
projection(x) <- "+proj=utm +zone=48 +datum=WGS84"
values(x) <- 1:ncell(x)


texts: Is it possible to add an ArcGIS attachment to a row in Oracle using oracle stored procedures?

texts:
So far I have been able to attach files to features within a featureclass using the python tools and ArcObjects.

Within Oracle I have been able to add features with geometries using the ST_Geometry functions from SQL Spatial however, I am wondering if there is a way to add attachments to these features within a stored procedure on the database as opposed to using python or ArcObjects


texts: Shapely's object.project method dealing with overlapping polyline

texts:
So I am interested in knowing how far along a polyline a given point sits. Python's shapely package has a method to do such a thing, see the below code:

from shapely.geometry import LineString

ip = LineString([(0, 0), (0, 1), (1, 1)]).interpolate(1.5)

print LineString([(0, 0), (0, 1), (1, 1)]).project(ip)
> 1.5

Now if I changed the polyline to fold back on itself it only returns 1 distance (the smallest distance), see the code below:

from shapely.geometry import LineString

ip = LineString([(0, 0), (0, 1), (0, 0)]).interpolate(1.5)

print LineString([(0, 0), (0, 1), (0, 0)]).project(ip)
> 0.5

In reality the answer to this question is 0.5 and 1.5.

My question is how do I get it to return all of the distances that may belong to the single given point?

UPDATE

I found that the shapely module uses PostGIS commands, specifically for the object.project method it uses the ST_InterpolatePoint command.


texts:
I calculated the distances between one set of xy coordinates and another set of xy coordinates.

So my near table has the ID variables, distance and rank.

However, when I join the near table to the csv file containing the xy coordinates, I lose all ranks that are not = 1. It only joins the closest distance.

How do I modify the joining process to include all distances and ranks?


texts: Rotating polygons by azimuth from attribute table using ArcPy

texts:
I want to rotate some polygons according to the value of azimuth from attribute table.
I try to use this code but there is an error that I can not understand.

CODE

import arcpy, math

fc = 'Polygons' #Change to match your data
rotatefield = 'azimuth' #Clockwise rotation in degress. Change name to match your data

#Function to rotate one input Point around a pivot point
def rotatepoint(point, pivotpoint, angle):
    #Source: https://stackoverflow.com/questions/34372480/rotate-point-about-another-point-in-degrees-python
    angle_rad = -math.radians(angle)
    ox, oy = pivotpoint.X, pivotpoint.Y
    px, py = point.X, point.Y
    qx = ox + math.cos(angle_rad) * (px - ox) - math.sin(angle_rad) * (py - oy)
    qy = oy + math.sin(angle_rad) * (px - ox) + math.cos(angle_rad) * (py - oy)
    return arcpy.Point(qx,qy)

#Rebuild each polygon with rotated vertices
with arcpy.da.UpdateCursor(fc,['SHAPE@',rotatefield]) as cursor:
    for row in cursor:
        polylist = []
        for part in row[0]:
            partlist = []
            for pnt in part:
                if pnt is not None: #Polygons with inner rings will have None pnt(s) which can be skipped
                    partlist.append(rotatepoint(pnt, row[0].centroid, row[1])) #Centroid is pivot point
                    polylist.append(partlist)
        row[0] = arcpy.Polygon(arcpy.Array(polylist))
        cursor.updateRow(row)

Error:

_Runtime error
Traceback (most recent call last):
  File "<string>", line 20, in <module>
TypeError: 'NoneType' object is not iterable_


texts: shapely wkt.loads - skip bad geometries?

texts:
I'm trying to convert an excel file into a Geopandas dataframe. I'm trying to use shapely and wkt to do so.

After reading the file into pandas, I have a geometry column. I'm trying to do the following to convert the geometry column to an 'object' data type:

from shapely import wkt

my_df['geometry'] = my_df['geometry'].apply(wkt.loads)

and I receive the following error:

WKTReadingError: Could not create geometry because of errors while reading input.

To try and diagnose what the problem is, I also tried to convert this excel file to a geodataframe, using the following code:

import geopandas

my_df = gpd.GeoDataFrame(my_df, geometry='geometry')

and I get the following error:

TypeError: Input must be valid geometry objects: MULTIPOLYGON (((1314112.145833299 1027703.927083313, 1314091.947916642 1027761.937499985, 1314232.583333299 1027811.447916642, 1314240.99999997 1027814.395833313, 1314246.739583299 1027794.468749985, 1314292.71874997 1027692.947916642, 1314282.18749997 1027689.010416642, 1314136.364583299 1027634.374999985, 1314112.145833299 1027703.927083313)))

It appears I only have one bad geometry? (Or it could just be the first bad geometry).

Can I just skip this bad one in the wkt.loads step?  I could not find documentation or other examples of how to do this.


texts: Why mNDWI can detect the shallow water area but ignores the deep water area?

texts:
I'm using Landsat 8 in Google Earth Engine to detect the water body in Jiangsu and Zhejiang area, China.

The formula is:

mNDWI = (Green-SWIR)/(Green+SWIR)

I applied this formula to GEE code, and it shows as follow:

As we can see, the blue area corresponds to the land, while the red area corresponds to the water body. The coastal line is clearly and correctly displayed. However, that only works for the shallow water area: As the distance to the coast increases, the red color fades away and the deep water area shows in blue.

The EVI result has the similar issue, I applied the formula:

EVI = 2.5*(NIR-RED)/(NIR+6*RED-7.5*BLUE+1)

And it shows as follow:

I promise that can't be cloud or some other disturb factors, I have already removed them before computing the indices. The shape of the abnormally displayed area is highly related to the shape of shallow water area, as shown below as a kind of yellow water:

I will also share my code here, hope anyone could figure out why it occurs, or is there any mistake I made, and how can I make the deep water area shown as what it should be.

//Choose Area
var region = ee.Geometry.Polygon({
  coords: [[[119.65, 36.28], [125.08, 30.63], [119.56, 29.65],[118.07, 33.81]]],
  geodesic: false
});

//Add image collection
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;

var non_mosaic = landsat8 .filterBounds(region) .filterDate('2017-01-01', '2017-12-31');

var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out cloudy shadow pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
  // Return an image masking out cloudy areas.
};

// A function to mask out cloud pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 5,5, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
};

//Will be used for frequency calculaton
var clouds_free = non_mosaic.map(maskClouds);

//display the mosaicked origional image and cloud free image
var mosaic_free = non_mosaic.map(maskClouds).median();
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};
Map.setCenter(121.07, 33.02, 6);
Map.addLayer(mosaic_free, visParams, 'Cloud free');

//The index results on a single day (Dec 31st 2017) are shown below

// Display the mNDWI.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an mNDWI image, define visualization parameters and display.
var mndwi = mosaic_free.normalizedDifference(['B3', 'B6']);
var mndwiViz = {min: 0, max: 1, palette: ['0000FF', 'FF0000']};
Map.addLayer(mndwi, mndwiViz, 'mNDWI');

// Display the EVI.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an EVI image, define visualization parameters and display.
var evi = mosaic_free.expression(
    '(NIR - RED)/(NIR+6*RED-7.5*BLUE+1)*2.5', {
      'BLUE': mosaic_free.select('B2'),
      'RED': mosaic_free.select('B4'),
      'NIR': mosaic_free.select('B5'),
});
var eviViz = {min: 0, max: 1, palette: ['0000FF', 'FF0000']};
Map.addLayer(evi, eviViz, 'EVI');


texts: Need to calculate field with special characters with ArcPy getting UnicodeEncodeError

texts:
I wrote a script in ArcPy and it attributes the features in a shapefile with the name of the municipality it's located in but some of the municipalities have special characters, Conceição do Araguaia for example, and causes this error - UnicodeEncodeError: 'ascii' codec can't encode characters in position 6-7: ordinal not in range(128)
My code looks like this:
for muni in muni_list:
arcpy.SelectLayerByAttribute_management(poly, "NEW_SELECTION", "{0} = {1}".format('Muni_ID', muni))
cursor = arcpy.SearchCursor(poly)
for row in cursor:
    print row.Muni_name
    arcpy.CalculateField_management(point, muni_field1, "{0}".format(row.Muni_name),"PYTHON_9.3","#")

The script runs fine until the last line of code is where the error occurs when it runs into "Conceição do Araguaia"


texts: How to make pgRouting understand bridges and tunnels?

texts:
I have I highway that runs on a bridge over a couple of blocks in a city. How do I get pgRouting to understand that cars can't move from the small blocks on the ground up to the highway on the bridge and vise verse?

What if I manually remove the nodes/vertices? Would a driving_distance understand that the highway and the smaller block-roads underneath aren't connected?

I'm using QGIS+pgRouting+Postgresql+PostGIS


texts: How can I prevent gdal_translate from rounding coordinates when clipping?

texts:
I want to clip a raster (a VRT of several Ascii files) to a set of coordinates. Ideally I'd like to use gdal_translate to do this, but I only get an accurate result when using gdalwarp. gdal_translate seems to round the clipping coordinates given in -projwin.

Here's a zoomed-in view of the corner of the clipped rasters when loaded into QGIS:

The corner of the clipped output should align with the diagonal line. The black rectangle (output of gdal_translate) is misaligned, the red rectangle (output of gdalwarp) is correct.

The gdal_translate command I used is:

gdal_translate -of VRT -projwin 429101.22 84985.077 429701.22 84785.077 83fb4946-fcdf-4536-b859-36b999efde7d_DTM.vrt one.vrt

QGIS shows that the extents of the clipped output are:

429100.0000000000000000,84786.0000000000000000 :
429700.0000000000000000,84986.0000000000000000

It looks as though gdal_translate has rounded or altered these.

The gdalwarp command I used is:

gdalwarp -dstnodata -9999 -cutline box.shp -crop_to_cutline -overwrite -of VRT 83fb4946-fcdf-4536-b859-36b999efde7d_DTM.vrt two.vrt

Where box.shp contains a rectangle with the same coordinates as those used in the gdal_translate command. QGIS shows that the extents of this clipped output are:

429101.2199999999720603,84785.0770000000047730 :
429701.2199999999720603,84985.0770000000047730

Which is much closer to the correct coordinates.

Is there a way to stop this rounding taking place when using gdal_translate?

Edit: I'm using GDAL 1.10.1

Edit: Here's the output of gdalinfo on a similar VRT:

Driver: VRT/Virtual Raster
Files: dd84ed63-b0c2-486f-a61c-8dea3d780552_DTM.vrt
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2883_DTM_2m.asc
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2884_DTM_2m.asc
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2885_DTM_2m.asc
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2983_DTM_2m.asc
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2984_DTM_2m.asc
   /Users/.../dd84ed63-b0c2-486f-a61c-8dea3d780552/sz2985_DTM_2m.asc
Size is 1000, 1500
Coordinate System is `'
Origin = (428000.000000000000000,86000.000000000000000)
Pixel Size = (2.000000000000000,-2.000000000000000)
Corner Coordinates:
Upper Left  (  428000.000,   86000.000)
Lower Left  (  428000.000,   83000.000)
Upper Right (  430000.000,   86000.000)
Lower Right (  430000.000,   83000.000)
Center      (  429000.000,   84500.000)
Band 1 Block=128x128 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999


texts:
I am processing a bunch of GPS data using a Qgis python script. The processing takes a LONG time to execute and when it is executing I notice that only a single CPU is pinned at 100%, the remaining CPUs are idle.

Before I get too deep into the details of my issue, I am trying to figure out how to create several processes to increase CPU utilization/processing speed for a task that is CPU-bound.

reading this question here: https://stackoverflow.com/questions/15143837/how-to-multi-thread-an-operation-within-a-loop-in-python I moved the contents of a loop in my script into a function then attempted to call it in using this form:

executor = concurrent.futures.ProcessPoolExecutor(10)
futures = [executor.submit(try_my_operation, item) for item in items]
concurrent.futures.wait(futures)

When I executed my script I got an error in the python error panel:

AttributeError: module 'sys' has no attribute 'argv'

After googling around a bit I found that using this at the top of my script should be a work around:

if not hasattr(sys, 'argv'):
    sys.argv  = ['']

When I execute my script now another instance of QGIS opens up and an error message with title "Private qgis.db" displays a message:

Creation of missing tbl_srs in the private qgis.db failed.
table tbl_srs already exists

Click OK and another error message witht he same title pops up:

Update of view in private qgis.db failed.
table vw_srs already exists

This then repeats its self 9 more times and I end up with 10 instances of QGIS open, no of which are actually doing anything...


texts: Annotation conversion into shapefile with appropriate Fonts Format in attribute table

texts:
My drawing file (.dwg file) is contained with Annotation with Gujarati Text in many type of fonts (e.g. Gujarati Saral, Nilkanth, Shruti...). After conversion of annotation to shapefile, attributes are seen in undefined characters.

So how can I get the same result in a shapefile?

Is there any process or solution to visualize the annotation characters in an attribute of shapefile.

[Tried to convert in ArcMap  from an attribute table > Properties > Appearance, which is temporary to see the characters]

Attachments:

1) Drawing File - Annotation Characters

2) shapefile - attribute table (Seen Undefined Characters)


texts: ol.control.MousePosition EPSG:4326 over 180/-180 longitude (OL4)

texts:
i will use a OpenLayers example in EPSG:4326, but the value of longitude is not wrapping on 180/-180 longitude.

How can i fix it?
My first problem is, i can't pick the values of latitude and longitude in a value for a if calculation.

if (lon > 180) {
lon = (((lon / 180) - Math.floor(lon / 180)) * 180);
}
if (lon < -180) {
lon = (((((lon * -1) / 180) - Math.floor((lon * -1) / 180)) * (180 * -1) + 180));
}

How do I extract the coordinates into a variable?
Or is there another solution?

The OpenLayers Example:
https://openlayers.org/en/latest/examples/mouse-position.html

The OL example in a JSFiddle:
https://jsfiddle.net/eqjjp4u8/


texts: Seeking tool to convert GTFS (General Transit Feed Specification) to SHP / KML?

texts:
I need to convert GTFS (General Transit Feed Specification https://developers.google.com/transit/gtfs/reference), data into some GIS data format like ESRI Shapefile or KML or something else.

Do you know any tool for this issue?


texts: Assigning values to lines around a reference point

texts: I am using ArcGIS 10.2 and I have a feature class with 150 lines. I want to assign values from 1 up to numbers of lines (150) to these lines based on the proximity of these points to a given reference point, the proximity should be based on the start or end point of the lines.


texts: QGIS 3.6 ordinary kriging produce one constant output

texts:
When I run kriging a+bx or any other model, I always get just one constant output. I have selected collect attribute.
Prediction and quality layer has just one value.

see output here:


texts: Changing projection for field calculator

texts: I've just succeeded in getting the area of a polygon using field calculator.  The polygon was drawn on a Google hybrid aerial photo. There are no units given in the calculation and I assume the area figure is based on latitude / longitude.  How do I go about changing the projection in order to get an area figure in hectares?


texts: Intersection of a vector linestring and underlying DEM

texts: Following is an screen-shot where ASTER DEM raster layer has been overlaid by a road linestring vector layer. I want to extract those DEM pixels which are lying exactly beneath the linestrings. I want to assign the corresponding elevation values to each and every linestring. Which PostGIS function should be used or any other approach? PostGIS ST_Intersection function for Raster Layers is not useful.


texts:
I need to create topographic maps for scales between 1:100,000 and 1:500,000. This is due to the fact that the official topographic maps in my country are offered in 5k, 10k, 25k, 50k, 100k and 500k. I really need to have a topographic data for cartographic work in 200k and 300k. I am looking at Maperative, Tilemill and OSM vector data, but would like to know if anyone has any experience or tips with this.

Thanks for any replies,

Robert


texts:
We are about to trial our GIS/LIS team moving onto AutoCAD LT as the functionality they use is minimal. I've installed it as a pilot but we seem unable to find how to connect to an SQL source.

Just wanted to confirm that it is even capable of connecting to an SQL source? Seems that the 'dbconnect' command is missing / does nothing?

AutoDesk has absolutely no mention of it on their site or support forum!


texts: Bounding Box approximation when calculating in euclidic space?

texts:
Although I will illustrate my problem using a database, I think this is a quite general GIS question. I hope I came to the right place to ask.

I am working with a system that runs on a normal relational database (without any SQL/MM functionality) and there is sadly no hope I can upgrade it anytime soon. But I need to maintain some basic geographic relationships regarding intersections of lines (Set A) and polygons (Set B).

Users may change individual lines and polygons and expect to get a list afterwards, regarding those that will somehow intersect. Now the exact intersection test isn't a problem and is dealt with in a Java backend. My current problem is that I have to check the cross product  of A and B anytime someone edits anything because the change could introduce a new intersection. No matter how small the change is, I can't tell in advance which exact relations could be affected.

So I would like to introduce some bounding boxes in the database, which hopefully will speed up the whole application because many obviously not related entities can be filtered out early.

I came up with the following structure and query to demonstrate the general idea:

CREATE TABLE bound_rect
(
  id serial NOT NULL,
  x1 double precision,
  y1 double precision,
  x2 double precision,
  y2 double precision,
  CONSTRAINT bound_rect_pkey PRIMARY KEY (id)
)

SELECT *
FROM bound_rect lhs, bound_rect rhs
WHERE lhs.x1 < rhs.x2
AND   lhs.x2 < rhs.x1
AND   lhs.y1 < rhs.y2
AND   lhs.y2 < rhs.y1

So this query is able to do some very general bounding box intersection tests (I am aware that things like the dateline are not taken into account, I can live with that). For me this seems like a similar calculation as it would be done in PostGIS' geometry mode. Performance is much improved, I am happy so far.

Now my actual question: Would I need to possibly "inflate" my bounding box if I want to make sure I don't get any "false negatives" when dealing with lat/lng coordinates on a sphere? I can live with "false positives", e.g. cases where the bounding boxes overlap but not the "real" entities. But if my query would filter out entities that do overlap because of too small bounding boxes, that would be a real problem.


texts: Convert a grid file from one set of coordinate systems to another

texts:
I have a grid file in a certain coordinate system. I'd like to change it so I have the same grid in geographic lat/long based on WGS84. This seems like something that gdal can tackle quite easily (maybe with gdalwarp or gdal_translate) but for some reason I can't find the proper way to do it.

I'm trying to run the following:

gdalwarp -s_srs '+proj=cass +lat_0=31.73409694444445 +lon_0=35.21208055555556
+x_0=170251.555 +y_0=1126867.909 +a=6378300.789 +b=6356566.435
+towgs84=-275.722,94.7824,340.894,-8.001,-4.42,-11.821,1 +units=m +no_defs'
-t_srs EPSG:4326 thick.grd thick2.grd

The proj4 definition taken from spatialreference.org. I get the following error message:

Warning: The target file has a 'grd' extension, which is normally used by the GSAG, GSBG, GS7BG, NWT_GRD drivers,
but the requested output driver is GTiff. Is it really what you want ?
ERROR 6: Unable to load PROJ.4 library (libproj.so), creation of
OGRCoordinateTransformation failed.

So I have two problems here. First, how do I let gdalwarp know what is the input format (in this case, Surfer grid) and the output format? Also, why can't it load the PROJ.4 library if I have it installed? Is there an alternative to PROJ.4 that will work on my system?


texts: Shapefile CRS damaged. How to bring it back?

texts:
I have two shapefiles with same location and extends (EPSG 4326). But one of the shapefiles CRS seems to be damaged after doing a lot wrong with reprojection... its not shown on the same position as the first, not damaged one.

Both have just to be on the same "place". How can I do that?

Damaged:

Good: