texts: How to convert a RGB tiff into pseudocolor with GDAL etc?

texts:
In QGIS I have made a singleband pseudocolor image out of one of the bands from a RGB GeoTiff. I choose one band and then make a linear color interpolation (see attached image).

I now want to do the same thing the command line way (GDAL/Mapnik/etc) so I can get a workflow going, but I don't really know where to start?

My idea is to split the RGB image into three greyscale files for each band and then give each band a pseudocolor map. Any ideas of how to do this?


texts: How to create a script using wgrib2 command

texts:
I have several model data with the following format

postwrf_d01_YYYYMMDD_INIT_f00000.gr2

where INIT = 0000, 0600, 1200, 1800, each consisting of up to several hours of forecast (f00000 to let's say f08400 with a six hourly interval).

Inside this gr2 data are multiple variables.

Now using wgrib2 utility, I can select just the variables I want using the following command

./wgrib2 postwrf_d01_20131107_0000_f00600.gr2 -match ':((U|V)GRD:10 m above ground|PRMSL|(ABSV|HGT|UGRD|VGRD|TMP):(850|700|500) mb):' -grib OUT.gr2

How can I create a bash script from the above wgrib2 command? I want the OUT.gr2 to have the following file formats.

1) d01_YYYYMMDD_ZZZZ_f0HHHH.gr2
2) YYYYMMDD_ZZZZ_f0HHHH.gr2
3) nwp.TC23.YYYYMMDDZZ.f0MMMM.gr2

For example
     input file : postwrf_d01_20160521_1800_f08400.gr2
    outfile shoud be
1) d01_20160521_1800_f08400.gr2
2) 20160521_1800_f08400.gr2
3) nwp.TC23.2016052118.f05040.gr2 ( 84 hr is 5040 min.)


texts: How can i add placemarkers (points) from click-event or database to openlayers3 map using geodjango

texts:
How can i plot place-markers in Geodjango using the methods:

a user click-event from the user on openlayers map preferably in reactjs
or alternatively from geodjango database that has PointField() attributes

Models.py:

from django.db import models

# Create your models here.
class Marker(models.Model):
    MarkerID = models.AutoField(primary_key=True)
    Name = models.CharField(max_length=40)
    mpoint = models.PointField()
    objects = models.GeoManager()
    Latitude = models.DecimalField(max_digits=10, decimal_places=6)
    Longitude = models.DecimalField(max_digits=10, decimal_places=6)

    def shpPoint(request):
        points = ShapefilePoint.objects.kml()
        return render_to_kml("placemarks.kml", {'places': points})

    def save(self,*args,**kwargs):
        self.Latitude = self.mpoint.y
        self.Longitude = self.mpoint.x
        super(Marker,self).save(*args,**kwargs)


texts: Regression Analysis with single value as predictor and grid as dependant variable

texts:
I'd like to perform a regression analysis in QGIS. Let's say for example that I have 3 grids as my dependent variable and 3 single values as my predictor. That is, for a single value of the predictor I'll have as many values as dependent variable as the number of cells in a grid. For another single value as predictor, I'll have as many values as dependent variable as the number of cells of another grid and so on. I want to perform a regression analysis between the predictor and the dependent, for example with a power law such as Y=a*(X^b) where X is a single value and Y all the values of a single grid. As a result, I'll have consequently two grids as output, one for a values and one for b values.
Is it possible to do so, maybe using Python?


texts:
I have written the following code to calculate the number of pixel classified as water within my region of interest.

I am manually adjusting the dates repeatedly to get a time series of flooded area/lake size i.e I am reading off the calculated "object" value and recording it in an excel file. My code has worked perfectly for 2019, 2018 and most of 2017 but fails for the years 2016 and 2015. I get this error:

Layer 1: Layer error: Image.lte: If one image has no bands, the other must also have no bands. Got 0 and 1.

Does anybody know how I fix this??

My code is as below:

//I am trying to construct a time series of flood extent over the duration of a whole year. The code will then be edited to be replicated for previous years as well

//This does all that I want.
//When you press run it will spit out two object values. The top one is the number of water pixles. The bottom value is the area classifiied as water
//Record these in excel spreadsheet.
//To get the next time step i need to go through and change the code each time i.e do this for each month, and then also for each year.
//So i need to change it all manually, rather than the code working it all out at for me

// Load Sentinel-1 C-band SAR Ground Range collection (log scale, VH, descending)
var collectionVH = ee.ImageCollection('COPERNICUS/S1_GRD')
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
.filterMetadata('resolution_meters', 'equals' , 10)
.filterBounds(roi)
.filter(ee.Filter.date('2016-01-01', '2016-12-31'))
.select('VH');
print(collectionVH, 'Collection VH');

//filter the date for the first half of january

var Jan1 = collectionVH.filterDate('2016-01-01', '2015-01-15').mean();

//smooth image to remove error
var SMOOTHING_RADIUS = 100;
var Jan1filtered = Jan1.focal_mean(SMOOTHING_RADIUS, 'circle', 'meters');

//mask image to just get water. This is using high certainity water mask
var Jan1mask = Jan1filtered.updateMask(Jan1filtered.lte(-23));

//map this layer

Map.addLayer(Jan1mask,{palette:"00FFFF"});

//The next step calculated the number of pixles classified as water (i.e they had a value of 1 rather than 0? Or 0 rather than 1, need to check)
// Sum the values of water pixels in the ROI. I did this one at a time and deletedthe old one after each step to save room i.e its been repeated over and over
var stats = Jan1mask.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e9
});
print(stats);


texts:
I'm trying to import an Excel sheet into QGIS3.

For that, I saved my spreadsheet as a .csv and in QGIS I went to "Layer->add Layer-> add delimited text layer..."

Unfortunately, I get the following screen.

Apparently, QGIS doesn't recognize the different fields from the original spreadsheet. Does anyone have an idea what the problem might be?


texts: Intersect of raster image on Google Earth Engine

texts:
I have this situation on Google Earth Engine: 3 different images, each of these are made by 17 different bands.
These bands represent the different year. I call these images with these letter: a,b and c. Image c is a binary map (1, no data). I would like to have a map that answer to this condition:

if b > b+1= a(b∩c) else a+1(b+1∩c)

for example if the values of a band of 2000 are major than 2001 the map respect this condition else l'altra.
I consider only the year and the next year (2000-2001,2002-2003 ....)

So I would if the condition is true the only pixel of a or a+1 image that respect the intersect condition.
I did some tests using only 3 bands (one per image) and I identified this possible solution in this command:

a.updateMask(b).updateMask(c)

But how I can do for respect my condition?

I think to an expression as the following code:

a.expression(‘?b>b+1=a.updateMask(b).update(c): a+1.updateMask(b).updateMask(c)’,
{‘a’:a.select(),‘b’:b.select(),‘c’: c.select() })

I obtain this error message:

Image (Error)
Image.parseExpression: Expression parse error at character 0:
'? b(y)>b(y+1) = a(y).updateMask(b(y)).updateMask(c(y)):  a(y+1).updateMask(b(y+1)).updateMask(c(y))'
 ^.

Where am I doing wrong?


texts: Alternative to ST_Union/ST_MemUnion for merging overlapping polygons using PostGIS?

texts:
New to PostGIS so I trust this is a very basic question but can't seem to find the right solution.

I have a materialized view consisting of ~4 million polygons. Some of them overlap. I need to merge those that overlap regardless of their attributes (the view only has two attributes; a unique ID and a geom column). Here's my query:

   CREATE MATERIALIZED VIEW dissolved_polygons AS
   SELECT ST_Union(c.geom) AS geom
   FROM overlapping_polygons c;

Every solution I've seen groups polygons by attribute, but when I remove the GROUP BY clause I get the following error:

  ERROR: array size exceeds the maximum allowed (1073741823) Query
  failed PostgreSQL said: array size exceeds the maximum allowed
  (1073741823) Here's my query:

After suspecting I was generating one massive multi-polygon I also attempted:

  CREATE MATERIALIZED VIEW dissolved_polygons AS
  SELECT (ST_Dump(a.geom)).geom FROM
  (SELECT ST_Union(c.geom) AS geom FROM overlapping_polygons c) AS a;

...but got the same error message. Another post here suggested ST_MemUnion might work but the table is way too big for such a slow process. Any ideas?

UPDATE 1

After testing with a smaller subset, I realize that the only way I can get the correct output is to use:

(ST_dump(ST_MemUnion(geom))).geom

...but I've executed this query on the real dataset and 9 hours later it's still running. I suspect it could take several days. Is there anyway to get away with using ST_Union without the memory overload, or is there a different command I could use I'm not aware of?

UPDATE 2

I've more clearly identified the problem but not the solution. According to (Much) Faster Unions in PostGIS 1.4 by Paul Ramsey, ST_Union (post PostGIS 1.4) aggregates all the data in memory before executing the union operation. The limitation of that memory storage appears to be 1GB, and my dataset is bigger than that. ST_MemUnion does what ST_Union used to do, which is union features row by row, without storing anything in memory. That solves the memory issue but it's incredibly inefficient and slow.
He talks about using ST_Accum() but I'm not clear on how/if this might help my query, given that I'm now using PostGIS 2.1.


texts: How can I find and fix broken raster metadata using QGIS?

texts:
I have a raster (TIFF) hillshade image that I would like to use in a mapping project. However, I'm having a hard time lining it up properly in QGIS. Through a series of investigations, I have determined that this problem might be caused by broken metadata. The raster was apparently created in ERDAS Imagine. I opened this raster in a text editor, to see if the metadata tags were visible. Sure enough, the following appears at the very end of the file:

#IMAGINE GeoTIFF Support
Copyright 1991 - 2005 by Leica Geosystems Geospatial Imaging, LLC. All Rights Reserved
@(#)$RCSfile: egtf.c $ IMAGINE 9.0 $Revision: 10.0 $ $Date: 2005/07/26 15:10:00 EST $
Projection Name = NAD_1983_HARN_StatePlane_Hawaii_3_FIPS_5103_Feet
Units = us_survey_feet
GeoTIFF Units = feet|IMAGINE GeoTIFF Support
Copyright 1991 - 2005 by Leica Geosystems Geospatial Imaging, LLC. All Rights Reserved
@(#)$RCSfile: egtf.c $ IMAGINE 9.0 $Revision: 10.0 $ $Date: 2005/07/26 15:10:00 EST $
State Plane Zone -5103
NAD = HARN

This data indicates the right projection, with the right units (feet) (this projection is functionally similar to EPSG:3760). ArcGIS is apparently able to read these tags, and readily assigns it the correct projection. However, when I view metadata on this file using QGIS or the QGIS browser I see that the projection of the raster is different. Most of the parameters are correct, but the units are set to meters. See screenshot from QGIS browser below.

Using ArcGIS, I removed the spatial reference of [a copy of] this raster, and then reset it using the Define Projection tool. Unfortunately, the metadata shown in QGIS remained the same. This leads me to two questions.

Where is this information coming from? Is there some other set of tags that QGIS is reading?
How can I remove or repair this broken metadata?

I'd prefer an open source solution, but I'll accept anything that works for now.


texts: Drawing from GeoJSON data using GeoServer and Leaflet?

texts:
I'm using geoserver and leaflet. I'm working to drawing from geojson data. But not happen. From geoserver data's coordinates are [668560.565133,4543124.14345],[668519.713326,4543116.052449] etc. I tried coordsToLatLngs, but not happen. Projection is EPSG:4326. With trying getFeatureInfo.

How can I convert or draw?


texts:
I tried to make Model in Model designer but I don't see my new fields in designer window, but only the offer of the existing layers.
What's wrong and what did I miss?
There is no connection lines.


texts: How to get rounded output from ST_AsLatLonText()?

texts:
I set a field with

UPDATE members SET latlong_degree = ST_AsLatLonText (the_geom)

But the result is (for example)

48°37’41.742”N 9°17’55.716”E

The seconds part is way too precise, what I want is

48°37’42”N 9°17’56”E

So I want to round the output…

Do you know a PostGIS command or other method to wrap the statement in, to accomplish this?


texts: Importing QGIS file in Eclipse for Android application?

texts:
Using ECLIPSE for android development, and also QGIS for customizing maps, I want to know if I can import the data in QGIS into my android application so that the data of my android applcation will be based on QGIS?

In other words, I want to have one database beneath both my QGIS and android applications.


texts: PostGIS raster table: How to get the global maximum (or minimum) value?

texts: Is there any straightforward and efficient way to get the global maximum (or minimum) value for a given band from a raster table in PostGIS?


texts: Using WMS getFeatureInfo() in shiny and leaflet

texts:
I like to get the getFeatureInfo by a click like it is mentioned in this post so that the information is shown in a popup window.
I tried trevis solution but it is not working for me. If I use addWMS() the WMS-layer is not visible and there is no info. If I use addWMSTiles() I don't get the getFeatureInfo Informations. Has anyone implemeted it correctly?
My code, after trevis post:
library(shiny)
library(leaflet)
library(leaflet.extras2)
library(XML)

ui <- fluidPage(

   # Application title
   titlePanel("Test"),

      mainPanel(
        leafletOutput("mymap")
      )
   )

server <- function(input, output) {
  output$mymap <- renderLeaflet({leaflet()  %>%
    addTiles(group = "base") %>%
    addWMS(baseUrl = "http://ows.terrestris.de/osm/service",
           layers = c("TOPO-WMS","OSM-Overlay-WMS"),
           # layers = c("OSM-Overlay-WMS"),
           group = "wmsgroup",
           options = leaflet::WMSTileOptions(
             transparent = TRUE,
             format = "image/png",
             info_format = "text/html",
             tiled = FALSE
           )) %>%
    addLayersControl(baseGroups = "base",
                     # overlayGroups = c("wmsgroup"))
                     overlayGroups = c("TOPO-WMS","OSM-Overlay-WMS"))
  })
}

# Run the application
shinyApp(ui = ui, server = server)


texts: How to add a shapefile programatically using PyQGIS

texts:
I am trying to add a Shapefile outside of the QGIS environment using PyQGIS. Ideally this would be done without creating a map in QGIS. I've started with the code below but am receiving the following error:

QObject::connect: Cannot connect <null>::raiseError< QString > to QgsVectorLayer::raiseError< QString >

Does the command below create a new map canvas?

from qgis.core import *
import qgis.utils

layer = QgsVectorLayer("F:\\IrrigatedLands\\FC_qgis\\boundary.shp", "testlayer_shp", "ogr")
#if not layer.isValid():
  #print "Layer failed to load!"


texts: What are Raster and Vector data in GIS and when to use?

texts:
What are raster and vector data in the GIS context?

In general terms what applications, processes, or analysis are each suited for? (and not suited for!)

Does anyone have some small, concise, effective pictures which convey and contrast these two fundamental data representations?


texts: Understanding GML structure for validation

texts:
I'm failing to convert this GML v3.2 XML and I'm not sure why:

<?xml version="1.0" encoding="UTF-8"?>
<gml:MultiGeometry xmlns:gml="http://www.opengis.net/gml" srsName="">
   <gml:geometryMember>
      <gml:Envelope srsDimension="2" srsName="EPSG:4326">
         <gml:lowerCorner>-67.4608695652174 -23.6347826086957</gml:lowerCorner>
         <gml:upperCorner>-32.0869565217391 -1.40869565217392</gml:upperCorner>
      </gml:Envelope>
   </gml:geometryMember>
</gml:MultiGeometry>

How is this supposed to be structured? I assume the envelope cannot be nested under the member tag.


texts:
I'm querying two different keys (to just make one API call) and would like to group them in the returned results. Is this possible or are the results always a per node list?

Code:

[out:json][timeout:25];
(
  node
  ["highway"="stop"]
  (around:5500,52.508088,13.465988);
);

(
   node
  ["highway"="traffic_signals"]
  (around:500,52.508088,13.465988);
);
out body;
>;
out skel qt;

The result should be:

  "elements": [
     [{stop},{stop},{stop}],
     [{traffic_signals},{traffic_signals},{traffic_signals}],
   ]

The current result is:

  "elements": [
     {stop},{stop},{stop},{traffic_signals},{traffic_signals},{traffic_signals},
   ]


texts: Set predefined scales in PyQGIS for standalone application

texts:
I am trying to figure out how to set the zoom scales to match those meant for OpenStreetMaps and Google Maps, etc., as I am using OSM for a basemap. I've looked all over the internet and can't seem to find any support for this particular goal.

I have an XML document (scales.xml) with all of the scales needed, but I just can't figure out how to implement it:

<qgsScales version="1.0">
    <scale value="1:591657528"/>
    <scale value="1:295828764"/>
    <scale value="1:147914382"/>
    <scale value="1:73957191"/>
    <scale value="1:36978595"/>
    <scale value="1:18489298"/>
    <scale value="1:9244649"/>
    <scale value="1:4622324"/>
    <scale value="1:2311162"/>
    <scale value="1:1155581"/>
    <scale value="1:577791"/>
    <scale value="1:288895"/>
    <scale value="1:144448"/>
    <scale value="1:72224"/>
    <scale value="1:36112"/>
    <scale value="1:18056"/>
    <scale value="1:9028"/>
    <scale value="1:4514"/>
    <scale value="1:2257"/>
    <scale value="1:1128"/>
    <scale value="1:564"/>
    <scale value="1:282"/>
    <scale value="1:141"/>
    <scale value="1:71"/>
</qgsScales>

I've tried reworking some of the code from thr QuickMapServices plugin, but I get no results from it:

def _load_scales_list(self):
    scales_filename = os.path.join(self.plugin_dir, 'scales.xml')
    scales_list = []

    xml_root = ET.parse(scales_filename).getroot()
    for scale_el in xml_root.findall('scale'):
        scales_list.append(scale_el.get('value'))
    return scales_list

def scales_list(self):
    if not self._scales_list:
        self._scales_list = self._load_scales_list()
    return self._scales_list

def set_nearest_scale(self):
    #get current scale
    curr_scale = self.mapCanvas().scale()
    #find nearest
    nearest_scale = sys.maxint
    for scale_str in self.scales_list:
        scale = scale_str.split(':')[1]
        scale_int = int(scale)
        if abs(scale_int-curr_scale) < abs(nearest_scale - curr_scale):
            nearest_scale = scale_int

    #set new scale
    if nearest_scale != sys.maxint:
        self.mapCanvas().zoomScale(nearest_scale)

This code neither throws an error nor does anything at all. Any ideas as to how I can get this working in a standalone window?


texts:
I add a listener to a button

 L.DomEvent
    .addListener(selectList, 'click', L.DomEvent.stopPropagation)
    .addListener(selectList, 'click', L.DomEvent.preventDefault)
    .addListener(selectList, 'click', this._onMarkerChange);

And then I want to move map by panto function.

_onMarkerChange: function (){
    this._handlingClick = true;
    for (var opt in this) {
        if(this[opt].selected)
        {
            var lat = this[opt].getAttribute( "data-lat" );
            var lon = this[opt].getAttribute( "data-lon" );
            this.L.Map.panTo( new L.LatLng( lat, lon ) );
            break;
        }
    }
     this._handlingClick = false;
},

It has an error --Uncaught ReferenceError: map is not defined
How should I do? pass map object into this event?


texts: Change the direction of the readers and writers in FME

texts:
Is it possible in FME to change the dataflow direction from:

Reading datasource A and writing to datasource B

to

Read datasource B and write to datasource A

And between the reading and writing you use many different kind of transformers. This meaning FME would have to do som kind of reverse engineer to be able to go the opposite way in all the transformers.


texts: Hyperlink (hotlink) creating multiple instance of QGIS

texts:
I have created a hotlink in QGIS 3.0, whereby the hotlink is designed to open another piece of GIS information.  The link works as I want, expect it opens a new instance of QGIS. How do i get the hotlink item to load into the existing instance of QGIS that is already running?

I have used the following options
Type =  Open url
Action [% "path"%]


texts: What is the algorithm used by ESRI for finding loops in a geometric network?

texts:
I've been trying to develop some code for ArcMap (VB .net) that would take a river network, identify and label the loops. These loops can be simple bifurcations or highly thread braided channels.

I was tinkering around with a sample of the NHD from the USA which is supplied as a geometric network. I tried the find loops Trace Task and was blown away by the fact that A) it did what I wanted and B) it is incredibly fast (less than a couple of seconds for a huge network).

Below is a screen shot of a bit of a network for Alaska with the polylines highlighted in red indicating they are part of a loop.

So I started looking on the internet to see if the algorithm is documented anywhere and I have failed to find anything.

I'm not a mathetician but I got the impression that I was looking at some sort of algorithm that detects  "strongly connected components in an acyclic directed graph". I have found many references to a Tarjans algorithm. I found this blog page which was the easiest attempt to explain these sorts of algorithms but as always they show a very simple graph which is not a realistic scenario when considering river networks. As the screen shot above shows you have loops within loops. I could not work out or even determine if it would work in a multi-threaded network.

So my question is what algorithm is ESRI using? What algorithm would find all the loops in a multi-threaded river network? I would like to try and implement it outside a geometric network as not everyone has an ArcInfo (Advanced) license and able to create a geometric network.

If anyone knows of a site\paper than describes such an algorithm in a simple manner that would be great!


texts: Adding Shapefile to Multiple Data Frames using ArcPy?

texts:
I'm trying to simply add a shapefile to 2 data frames of a Map Document, the shapefile is already in a third dataframe of the same Map Document. I was using the guide below but am getting the following error: ValueError: Object: CreateObject Layer invalid data source

Adding shapefile or feature class as layer in ArcGIS Desktop using Python/ArcPy?

import arcpy
from arcpy import env
import arcpy.mapping
env.workspace = "C:/Users/Chris/Documents/501_Customization                               /lab4/E_ercise_10_export/Austin"

#Get Shapefile
parkslayer = arcpy.mapping.Layer("C:/Users/Chris/Documents/501_Customization/lab4/E_ercise_10_export/Austin/parks.shp")

#Get Map Document
mxd = arcpy.mapping.MapDocument("CURRENT")

#Get data frames
df = arcpy.mapping.ListDataFrames(mxd, "*")[0, 1]

# add the layer to the map at the bottom of the TOC in data frame 0 and 1
arcpy.mapping.AddLayer(df, parkslayer, "BOTTOM")

# Refresh things
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd, df1, parks

print "Shapefile successfully added to data frames"


texts: Spatial join by argument closest, but based on majority of vertices

texts:
I'm often working with different polyline datasets which represent street networks. One of the most common tasks I have is to join information from one layer to another. Often there are no mutual attributes to utilize.

Geometries between layers are similar, but they're far from identical. First thought would be to use spatial join based on argument closest. Alongside the fact that the argument doesn't exist in core QGIS function, as far as I know it is based on a single, closest vertex. Since we're talking about wide road networks, the single, closest vertex might not be the correct one (see image). It's more about two lines sharing approximately the same path.

In the image, based on approximate paths of the three lines shown, the orange line should be joined with the blue one. Yet, since a vertex of the green line is the closest to any of the vertices of the orange line, the join would actually happen between them.

So I reckon the solution should take into account all of the vertices, get their closest counterpart and do a join based on their majority.

My question is: are there any existing open source tools or combination of them that could tackle this efficiently or should I start planning a script?

I'm using QGIS 3.4 on Windows.


texts:
I've got a dxf in AutoCAD with a different Z value in each vertex. I'm importing this dxf to QGIS and running the LRS plugin to work with some chainages.

My problem is that when running the chainage I'm not getting the same than my peers are getting in ArcGIS, by as much as 15m.

I think the issue could be that LRS is not reading the 3D elevations and treating the line as a 2D. Although if this was the case, the difference should be quite highers, since it's a 450km and there is more than 300m of difference between elevations at start and end.

Anyways, how do I make sure that the LRS plugin is taking in consideration the Z value?


texts: Creating points along line at ellipsoidal distance

texts:
I need to create points along a polyline in QGIS 3 at every 1000m of ellipsoidal distance and have tried several tools, as suggested in other answers:

Points along geometry
QChainage plugin
Convert lines to points (SAGA)

All of these can create points along a straight, two-node line, but do so at planimetric distances, so when I measure with my measuring tool (in planimetric mode) I get 1000.000 m:

However, if I switch on the WGS84 ellipsoid distance measurement, I get a deviation:
I'm working in WGS84/UTM32N. I already know the ellipsoidal length of my line and the last measurement point at the end of the line should have the same ellipsoidal value.
Is there any way to do this in QGIS 3? Am I missing some basic coordinate system logic? I am working in a UTM CRS, which is planimetric, I guess?


texts: Vector to raster conversion: Gaps within the vector boundary

texts: Is there a tool in ArcGIS that will make the output raster completely encompass the vector boundary? Somewhat similar to select by attributes where all overlapping features however minuscule is added in the selection. What I do now is create a buffer around the polygon with the distance same as the cell size of the output raster I want, but this is an extra step which if has a proper tool (or setting) counterpart will greatly streamline my workflow.


texts:
US census summary file information is available on the census ftp server: ftp2.census.gov/acs2013_5yr but I can't seem to find the data profiles information anywhere.

The only source appears to be the census data api (census.gov and api.census.gov).

Is there any place you can download this raw data or is this API the only source? I've scoured the census site and I hope I'm not missing something obvious.


texts: custom equal area projections with proj4

texts:
I need to define equal area projections for various countries around the world (in proj4 terms). I’ve taken a stab and designing them myself but would like to verify that my method is correct.

For example, for the part of the world encompassing Borneo, Malaysia and Papa New Guinea, I’m modifying the ESRI: 102028 projection for South Asia using the 1/6 rule.

The lat span of interest is 6.7 to -11.82 degrees so a range of 18.52 deg.
1/6 of this is 3.08.

lat_1 = 1/6 south of my northern extreme = 6.7 - 3.08 = 3.62

lat_2 = 1/6 north of my southern extreme = -11.82 + 3.08 = -8.74

I would also readjust the central longitude to 117 degs.

So rounding a bit, my new projection is:

"+proj=aea +lat_1=4 +lat_2=-8 +lat_0=0 +lon_0=117 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

Can someone confirm that I’m thinking about this correctly?


texts: Showing description of value map in print layout in QGIS 3

texts:
How can I show the description of a value map (joined layer) in a print layout in QGIS (QGIS 3.4.15 or 3.10.6). It just shows the value...

Data is Geopackage, point layer joined with table (no geometry). I select the values in the form with a value map (value and description, both with different entries). In the form and the attribute table QGIS shows the description as it should be. In the print layout (using atlas) it shows the description as well, but if I use a joined field it just shows the value, not the description.


texts: Converting vector points to geo coordinate points

texts:
So I have an SVG file that has paths for a floor in a building.

I want to convert these polygons into something like GeoJSON so that I can draw them in Mapbox.

I already know the geo coordinates and dimensions of the building Im trying to map.

To start I've been following the solution in This Answer

I got as far as creating the shapefile image but I'm stuck after that.
The Affine transformations plugin is incompatible with my version of QGIS3.

Can you suggest another way I can try to convert my SVG polygons into Geo Coordinates?


texts: How to extend QGIS Graphical Modeler with Python Scripts

texts:
I'm trying to extract lat longs from my centroids, then feed the lat longs to the Google distance matrix API using another python script. I exported the graphical model as a python script:

##Drive Time=name
##GridSpacingInMeters=number1000
##Extent=extent
##CentroidsWGS84=output vector
##GridWGS84=output vector
outputs_QGISVECTORGRID_1=processing.runalg('qgis:vectorgrid', Extent,GridSpacingInMeters,GridSpacingInMeters,0,None)
outputs_QGISREPROJECTLAYER_1=processing.runalg('qgis:reprojectlayer', outputs_QGISVECTORGRID_1['OUTPUT'],'EPSG:4326',GridWGS84)
outputs_QGISPOLYGONCENTROIDS_1=processing.runalg('qgis:polygoncentroids', outputs_QGISREPROJECTLAYER_1['OUTPUT'],None)
outputs_QGISEXPORTADDGEOMETRYCOLUMNS_1=processing.runalg('qgis:exportaddgeometrycolumns', outputs_QGISPOLYGONCENTROIDS_1['OUTPUT_LAYER'],0,CentroidsWGS84)

Then I added more lines to the above to create a list of the coordinates, but I'm getting a IndexError: list index out of range

name = 'CentroidsWGS84'
layer = QgsMapLayerRegistry.instance().mapLayersByName(name)[0]
features = layer.getFeatures()
attributes = []
coords = []

for feature in features:
    atributes.append(feature.attributes())

for attribute in attributes:
    coords.append(attribute[5:7])

I suspect this is because the layers that I'm creating on the graphical modeler is not created yet. How do I solve this?


texts: Using quotation marks in ModelBuilder of ArcMap

texts:
The syntax in the ModelBuilder of ArcMap varies when trying to use the 'Name' output.
When do I use double quotes, single quotes, or no quotes?


texts: Creating shapefile based on map with ArcGIS Desktop?

texts:
I'd like to create a shapefile with ArcGIS of the map attached below. Published in an econ paper, the map divides the territory of the German Union (post 1815) in 19 states (more precision can be found here https://economics.mit.edu/files/6330)

Quite new to GIS, I am looking for a strategy that could allow me to create a shapefile of the German Union with identical borders/internal state delimitations.

Harvard's HGIS department produced some shapefiles at Germany States Boundaries.

Thus I thought that a good strategy would be using Harvard's purple shapefile with States Boundaries to obtain the desired shapefile, with the desired partition of 19 states.

Unfortunately, turns out that the desired partition of 19 states is finer than the one in Harvard's shapefile. I would need to draw additional internal borders / divide some polygons into multiple ones (ex East) -- and I don't know how to this, despite having searched on the web for hours.


texts: Installing updates/patches for bug fixes in QGIS?

texts:
I'm trying to fix this bug in QGIS 2.0.

I suppose the code, under the "Files Changed' tab needs to be copied/installed somewhere?   Could you please advise on how to do so? I cannot find any help page on this - but may have been Google-ing the wrong terms.

More broadly (but not my immediate concern) - why are these changes just not corrected in the original source code?


texts: Is there a way to download and edit an OSM file in QGIS?

texts: I am working with OSM maps in QGIS.  I am looking for a way to extract a shapefile of roads, or buildings, from OSM so I will be able to edit them.


texts: Is there an exhaustive, searchable catalog of all GIS web services (e.g. OGC WMS/WFS, REST, SOAP)?

texts:
Many of us have struggled to get a list of GIS services that are available on the internet. Specifically map services that expose OGC based interfaces like WMS, WFS or Web Service interfaces like REST, SOAP etc.

For instance, just for oceanic and marine data there are several online catalogs such as
1) FDGC
2) NOAA NowCoast
3) IMO Mest
4) NASA

If you start with weather, transportation, demographics, agriculture, this list becomes huge!

These sites are not easy to use, don't have the most updated info and don't have a good taxonomy/ classification of services towards easy discovery of such services.

Is there a better attempt out there to create an exhaustive, searchable catalog of all such services? I would like to know their usage rights, free or pay etc. Uptime, performance etc. would also be nice to know.


texts: How to transform .qml/.sld in .xml with QGIS

texts:
I have a bunch of .qml and .sld files. I want to convert them to .xml, but I still havent found a proper way to do it.
Does anyone know any plugin or website that I can do this conversion?


texts: QGIS Rastercalculator: unexpected results for cos() function

texts: I've calculated the slope values of a DEM and created a raster with values from 0 to 51,3 (degree). I wanted to calculate the Cos of these pixel values with the Raster calculator (to calculate in the next step the total area) expecting only positive values (cos(0...90) should be positive) but around 50 percent of the values turned out negative and I can't figure out why.


texts: Does reduceRegion take care of sticking the neighboring tiles?

texts:
Let's say we pull an image on google earth engine where the entire region is spread across two neighboring tiles of Sentinel. i.e. half of the region is in one tile and half of it is on the other tile.

When we do reduceRegion, does Google Earth Engine automatically take care of such event? Or do we need to mosaic/stitch the tiles first and then do reduceRegion?


texts:
I have a feature class that contain many overlapping polygons. Each polygon has a field value which I would like to add together wherever the polygons overlap. There are about 100 polygons each with complex boundaries and some are multi-part polygons. Below is conceptual diagram of what I would like to do.

I know I can do this by converting each polygon to raster and adding them. The problem is that since their boundaries are complex, it takes a very long time to convert. I was wondering if there is a solution without converting to raster?


texts: reproject vector layer with pyqgis with new CRS does not worl

texts:
In QGIS I am running a script in th python toolbar to export all visible vector layers to a GEOJSON with this snippet:

canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()
for i in allLayers:

    if i.type() != 0 :
        print(i.name() + " skipped as it is not a vector layer")
    if i.type() == 0 :
        exp_crs = QgsCoordinateReferenceSystem(4326,QgsCoordinateReferenceSystem.PostgisCrsId)
        print i.crs()
        print exp_crs
        qgis.core.QgsVectorFileWriter.writeAsVectorFormat(i, '/home/' + str(i.name()) + '.js', 'utf-8', exp_crs, 'GeoJson')

If my vector file has a different CRS it keeps using the vector file CRS and doesn't use the defined CRS EPSG:4326... What am I doing wrong?
The whole project context can be found on GITHUB


texts: Problem importing qgis.core inside docker container

texts:
I'm trying to use QGIS libraries in docker container. When I try import qgis.core I get this error, really I can't figure out what is the cause.

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/qgis/__init__.py", line 36, in <module>
    from qgis.core import QgsFeature, QgsGeometry
  File "/usr/lib/python2.7/dist-packages/qgis/core/__init__.py", line 34, in <module>
    from qgis._core import *
ImportError: No module named _core

I've added proper paths to PATH environment.

sys.path.append("/usr/share/qgis/python/plugins")
sys.path.append("/usr/share/qgis/python")
sys.path.append("/usr/lib/python2.7/dist-packages")

and this is my PATH

/usr/share/qgis/python/plugins:/usr/share/qgis/python:/usr/lib/python2.7/dist-packages:/usr/local/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:


texts: PCS and GCS aligned but XY coordinate values not. ArcMap .dwg to .shp conversion

texts:
I have converted a .dwg to a .shp using the CAD2SHP app. I made sure that there would not be a .prj file when I did the conversion. Defined projection once I opened up the file in Arc as a .shp. Then I bought in another .shp, with the same GCS and PCS. I don't see any differences between the two, and yet they will not line up.

Arc has not shown me any errors in the process of  converting or projecting any of the data.

The only clue I can provide is that they have different coordinate values, one in the thousands and the other in the millions. In full extent they appear quite far apart.

Any insight - or a link to someone with a similar issue.


texts: Pause ArcPy script until feature class lock goes away

texts:
My script needs to copy data into a feature class to an SDE geodatabase that might have a lock on it from another application. I shouldn't break the lock (that app has priority over a data update), but I also don't want the script to quit just because it tried at the wrong moment.

How can I set it to keep trying until there is no lock? I'm looking for something like a try/except, but more like try/keep-trying-every-15-seconds.

I considered making a separate function and adding a pause:

try:
    funcCopyData()
except:
    pause 10 seconds
    funcCopyData()

but wouldn't this still be limit me to only two attempts?


texts:
I have two shapefiles, one with points and one with polygons.

I have written a query for checking whether the point is inside a polygon and if it is inside, assign a value.

However, I want to assign the value 0 if the point is not inside the polygon.

Here is the query:

SELECT polygon.gid, time, avg(value)
FROM   polygon,geo
WHERE  ST_Contains(polygon.geom,geo.geom)
GROUP  BY polygon.gid,time

I have multiple values for each point at different times.


texts:
I am using OpenLayers to create a map and I am using:

map = new OpenLayers.Map('map')

to create the map layer. But the map appears somewhere down in the web page. How can I tell where the map is to be placed on the page. I read the OpenLayers documentation and found that there are some attributes that can be defined with the above map variable like:

map = new OpenLayers.Map('map', {
    maxResolution: 2
    });

But I am not able to see any difference in position of the map. Any suggestions?


texts: QGIS running processing algorithm in Python Console does not store Outputfile?

texts:
I have a process which I wish to automize with a Python Script. This Python Script accepts 2 input layers (1 point layer and 1 line layer; both layers do not have the same CRS). The script should then calculate the nearest line to each point.

My solution for this includes 3 steps:

1st step: reproject input layer (the line layer to be exact) with the processing algorithm "qgis:reprojectlayer" (in order to have two inputlayers with the same CRS)

2nd step: Calculation of the nearest lines to each point with the NNJoin plugin

3rd step: Transfer attribut with the field calculator

Before I started writing the Python script, I opened the Python console to test the commands I would need. First of all, I tried to test the qgis:reprojectlayer algorithm (for that, I did a "import processing").
And then, I checked the processing history and copied the usage of the algorithm into the console (I have run the algorithm before in order to do that).

And here is the problem: If I run this algorithm from the toolbox, it generates a output file which is instantly loaded into qgis. The same happens if I run the algorithm directly in the history window (with a doubleclick). But whenever I copy the functioncall (from the history) into my python console and run the command, it looks like this:

It says OUTPUT: memory:, but it does not load anything into qgis and I have also tried it with a different output-file-path but the output folder remains empty. It looks like the algorithm does not generate a file?
And how can I assign the output file (which should be generated by the reprojectlayer-algorithm) to a variable within my python script?

My Python Script looks like this at the moment (I defined the inputs and output):

Any suggestions?


texts: Cannot install plugin due to Bad Gateway error

texts:
While attempting to install a plugin via the Manage and Install Plugins window I receive the following error:

> Plugin installation failed:Error downloading
> https://plugins.qgis.org/plugins/quick_map_services/version/0.19.10.1/download/
> - server replied: Bad Gateway

Also, when I attempt to follow part or the full link in this error I receive an message that reads 502 Bad Gateway nginx/1.4.7in my web browser. I also tried installing a different plugin and received the same error. I just upgraded from QGIS 2.10 to 2.18. Am I doing something wrong?


texts: Anybody know how to add an arrow style to a EE LineString?

texts:
I am trying to add an arrow style to a LineString in Earth Engine, pointing in the direction that the LineString coordinates were created.
The only style options I see are solid, dotted, dashed.  Does anybody know how to do this?


texts:
When I try to transform geometries from SRID 2321 to 4326, the transformed coordinates are swapping. Why is this happening?

Below are the details

Source

SRID - EPSG:2321

Source Data Projection Definition:

WGS84 Bounds: 31.5000, 35.9000, 34.5000, 42.1000
Projected Bounds:364573.0642, 3975550.1410, 635426.9358, 4663916.1654
Area: Turkey - 31.5°E to 34.5°E

below are the sample points :

470343.6,4402310.3
470105.0,4402410.8
470236.5,4402700.1

Destination

SRID - EPSG:4326

Below are the sample points after conversion:

32.6539666652,39.753664831
32.6511781435,39.7545616174
32.6526994043,39.7571716999

When I try to project these point on a Google map they are referencing to Iraq when they should refer to Turkey.

But when I tried with reversing the coordinates of the above points then they are referencing to correct locations on Google maps

Why is this happening? Do I miss anything?

EDIT:

Software used for conversion : Spatialite v 4.1.1 and this conversion is
confirmed using this link.

Actually my source data is from ESRI personal geodatabase .mdb file. I used QGIS to import them into Spatialite and when I try to project them on google/Apple maps all data is referencing to a different location.


texts:
With QGIS 3.6.2 on Ubuntu 18.04 I can't get the startup.py file to execute. Checked the docs and default folders with QStandardPaths.standardLocations(QStandardPaths.AppDataLocation) and added the startup into the final folder /home/jaume/.local/share/QGIS/QGIS3/profiles/default/python nothing happens.
The test code I want to execute is:

f = open("/home/jaume/aaa.txt","w+")
f.write("This is a line")
f.close()

I've tryed other silly code, added execution permissions, tested that this code runs in the python3 console, etc. Nothing happens. Any advice?


texts: Selecting countries that are being crossed by a river

texts:
Can you please advise how to mark countries that are being crossed by a river?
Meaning, to mark it not manually, but with a sort of function that I can later export this layer? I tried using Select by location... with no success. Created a buffer around the river. Whatever is in that buffer should be selected as the crossing country.
I am using QGIS.


texts: ArcGIS Pro - Open .aprx files to edit file paths

texts:
Two times a year a positional check of railway assets is performed in ArcGIS Pro for 21 different areas by receiving updated data. Every half a year in each area new data is automatically prepared in a standardized structure and needs to be loaded in a new ArcGIS project file (.aprx). Therefore the file paths in the .aprx file should change to the new data based on the area and date. I know the .aprx files can be unzipped to show .xml files, but I am not sure of best way to edit all the filepaths to load the new data.

This problem is very similar to editing file paths in QGIS .qgz files.

How can I change the file paths in the .aprx file easily?


texts:
I have a few CH1903/LV03 layers (lake, river and the country border of Switzerland) and a layer in WGS84 (two points of two cities, Vaduz and Bregenz).

On-the-fly-transformation is activated but Vaduz (Liechtenstein) and Bregenz (Austria) are still shown in Switzerland. Data source is a PostGIS database. Any idea why?


texts: I am looking at switching over a computer from Windows to Ubuntu.  I have been a Windows user all my life and been using GIS for over a decade.  I am curious about the migration to Ubuntu.  Mostly, what packages, bindings etc I will need/want to get.  With Windows I have downloaded Python and the bindings, GDAL, FWTools etc etc.  Will it be the same with Ubuntu? Or are packages bundled differently in Ubuntu?


texts:
I'm trying to create a topology of the public transport in Paris. I use a OSM export, osm2pgrouting with the correspondent mapconfig.xml and Postgis.
Everything works fine, except one thing.

In Paris, there are two kinds of public transportation on railways: the metro and the "RER", which is basically a metro that go further in the suburbs.

I'm am able to import both networks (metro and RER), but I'm not able to connect them properly, i think because on OSM the shared stations are mapped as relations of two nodes (one station for the metro and another one for the RER while it is the same).
One example can be seen there:
http://www.openstreetmap.org/relation/2154932

According to :
http://workshop.pgrouting.org/chapters/osm2pgrouting.html
It seems that osm2pgrouting doesn't handle relations.

Is it true that osm2pgrouting doesn't handle relations, or am i over-interpreting ?

If this is true, how can i manage to overcome this problem ?


texts: Google maps white striations / false artifacts

texts:
What are the common white striations running from NW to SE visible when you zoom in around 5mi scale? For example at (+81, -80)? But they appear 'all over'. Are they snow drifts? or satellite sweeps?

Generally: Where can I find documentation on Google Maps false artifacts?


texts: How do I .replace/remove within labelling

texts:
I would like to remove N_A and N123#N123 from my Labels. How do I do that?

So far:

str([RTN]).replace('N_A',' ')

I still Need to remove the eg. N123#N123 as well from the Labels, but keeping the single eg. N123 Labels.

Attribut Tabel:
N_A        to be removed
N123#N123  to be removed
N123       must not be removed
N321       must not be removed


texts: Adding click event to every polygon on map in leaflet to use map.fitBounds()

texts:
Here is my code. But it's showing only last polygon which is iterated from array. (For binding polygon I am using array)

var polygon, i;
for (i = 0; i <= latlngs.length; i++)
{
    if (i < 2) {
        polygon = L.polygon(latlngs[i], { color: 'green' }).addTo(mymap);
    }
    else {
        polygon = L.polygon(latlngs[i], { color: 'red' }).addTo(mymap);
    }
    polygon.on('click', function (e) {
     mymap.fitBounds(polygon.getBounds());
    });
}


texts: How can I improve my expression including an aggregate-function?

texts:
I have a polygon layer and a point layer ("wka"). I want to count the points in each polygon. I know I could do it with the following tools: "count points in polygon" or "join attributes by location". This works, but the disadvantage is the geodata management I have to do: save in a new location, import in the project, copy the styles form the old polygon layer to the new, etc.
I found out I can use the aggregate function in the field calculator. I tried that expression and it surprisingly works :)
aggregate(layer:= 'wka', aggregate:='count',
expression:="fid", filter:=contains(geometry(@parent), $geometry))

I am wondering is there a "cleaner" way to write that expression? Does the part expression:="fid", make any sense? I could write any attribute name for "fid" and it still works. But when I delete this part from the expression, the expression doesn't work anymore. I am wondering why.


texts: scripted new items don't show up in the "items" panel on Print Composer

texts: I'm working on a plugin to put items (pictures, arrows, labels) on the Print Composer layout. Everything is fine with the exception that once they are drawn, none of these show up in the "items" panel until saving and reopening the project, or duplicating the composer (only in the duplicate). It is not a terrible problem but it is not ideal either. Does anyone know why and how to solve this?


texts: How to resolve the error while downloading MOD10A2 image using MODIS package in R?

texts:
I am trying to download the MOD10A2 image using MODIS package in R
However, each time it is failing with following error

Error in sf::gdal_utils(util = "warp", source = gdalSDS, destination =
ofile,  :    gdal_utils warp: an error occured

I am trying this
library('MODIS')
library('gdalUtils')
library('sf')
EarthdataLogin(usr = USERNAME, pwd = PASSWORD)

MODISoptions()
lap = 'I:/'
odp = file.path(lap, "MODIS")
MODISoptions(localArcPath = lap, outDirPath = odp)

date = transDate('2000-02-24')

runGdal(
  'MOD10A2',
  begin = date$beginDOY,
  end = date$endDOY,
  tileH = 24,
  tileV = 5,
  job = 'SnowCover',
  SDSstring = NULL,
  overwrite = FALSE,
  maskValue = NULL,
  outProj="GEO"
)


texts:
I constructed a python query to look for all the bus routes passing by a given box. However, I need to speed up the query as much as possible. I only need the lat/lon coordinates of each node representing the bus stops and the name of the routes that they belong to. My current code:
import requests

overpass_url = "http://overpass-api.de/api/interpreter"

bbox = [48.87542724909715, 2.1707384550740683, 48.88884184835508, 2.1821696229817267]
bbox_str = '('+str(bbox)[1:-1]+')'

overpass_query = '''
[out:json];
(
relation'''+bbox_str+'''['route'='bus'];
);
out body;
>;
out skel;
'''
response = requests.get(overpass_url, params={'data': overpass_query})
data = response.json()


texts: Algorithm to compile Lisiecki stack? (5MM-yr Temperature Record)

texts:
I have obtained the Lisiecki stack data from www.lorraine-lisiecki.com/stack.html

How do I convert the d18O measurement into a figure for average surface temperature? What is the best algorithm?

I found this on wikipedia: If the signal can be attributed to temperature change alone, with the effects of salinity and ice volume change ignored, a δ18O increase of 0.22% is equivalent to a 1 °C (1.8 °F) cooling.

Any help would be useful. I am also trying to find the dataset for the fabled 65-million-year temperature record.


texts: Calculating roundness/compactness of polygon?

texts:
I am trying to find ways to quantitatively describe the shapes of various polygons. For my project, these polygons represent lakes, rivers, lagoons, and parks. So they can be almost any shape. One easy metric is to calculate perimeter vs. area, which is at best an only slightly useful metric. But I would also very much like to be able to say something about the 'roundness' of a polygon. Or how 'compact' the shape is on a map.

The only way I can think of doing this easily is to calculate the area of each polygon in relation to a bounding box for that polygon (which I already have). But this seems like a poor solution.

So now I am thinking of something more like this - take the centroid of the polygon, add on a series of buffers of increasing areas (say 50%, 100%, 150%), then compare how much overlap there is between each buffer and the original polygon. A perfect circle will have perfect overlap at 100%, and I can use the 50% and 150% buffers to judge how much and in what way each polygon differs.

But even that feels cumbersome, and like a poor workaround for what somebody else has probably already figured out far better.

For reference, at a minimum I will need to be able to look at the resulting indices for the shape of various polygons, and be able to make an educated guess as to their source (River? Reservoir with dendritic shape? Lake/Lagoon? Park?)


texts: Merging features from different layers by field value in QGIS

texts:
I have various layers with polygons which each represent a land use class (e.g. "roads", "vegetated terrain", "water", etc.) These layers have fields which further classify the features' land use (e.g. "bike path", "forest", "lake").

I'd like to re-classify information from all layers into new categories ("paved area", "roofs", and others). I will be combining features with different field values from various layers for that.

Is there an easy way to do this in QGIS? So basically:

Layer 1, field subclass = X
Layer 1, field subclass = Y
Layer 2, field subclass = Z

should be merged into a new shapefile representing the output class I want.


texts:
I have a 23-band .sav file that I need to load in ENVI.  The 23 bands are 23 time points of MODIS NDVI data.  I ultimately want to run a classification on this data.  I have been using Data Manager to load the .sav file, but have only been able to load three bands at a time (RGB).

How can I load the ENTIRE 23-band file at once?


texts: Rule-based classification not working on QGIS

texts: I'm following the official QGIS tutorial on the QGIS website and I'm trying to create a rule-based classification as in the tutorial and I've followed the same steps as they have but it's simply not working. The whole polygons in the file I'm trying to classify only changes color to that of the default (no filter) criterion. I've attached some pictures here. The red arrow on the first picture points to the default criterion I mentioned earlier. What am I doing wrong?


texts: Public shp file showing Berlin subway entrances?

texts:
I was inspired by someone else's work to try finding the Berlin property farthest from a subway station. In order to do that, I will need two things:

First, I need a shapefile of all the building addresses in Berlin. I believe I have found this here:
http://opendata.esri-de.opendata.arcgis.com/datasets/d4ac5d6ff99140819d2033038664096c_0

I will also need a shapefile showing all the subway (U-Bahn) entrances in Berlin. Does anyone know of such a shapefile?

Bonus heaps of gratitude will come your way if you can dig up a shapefile of S-Bahn stations as well, though this is not strictly necessary. Note that I am looking for the entrances/stations, rather than the lines themselves.

Assuming I can locate a shapefile showing U-Bahn lines, I will be doing my analysis using python and QGIS. Here is the original inspiration for this project:

http://iquantny.tumblr.com/post/99544282749/found-the-manhattan-apartment-thats-the-farthest


texts: Positioning coordinates in middle of map data tile?

texts:
Above you can see an image tile which is generated based on a lat, lon position. As you can see that the generated image and data tile contains a lot of unnecessary mountain region. For getting the tile image, please refer to the following link:

https://developer.here.com/rest-apis/documentation/enterprise-map-tile/topics/key-concepts.html

Is there a way to generate a tile with the provided lat lon being in the middle of the tile?


texts:
I have a shapefile that contains elevation data.  I want to eventually convert this to a slope layer.  I'm assuming I have to convert to a raster, then a DEM, then slope.

Is that correct?

If so, could someone point me in the right direction?

I'm not very familiar with the program.

When I try to use the inverse distance weighted method, I get a layer that has one point only. It cover the entire area of my original layer and more.


texts: Help in field calculator expression editor hidden

texts: I don´t have in field calculator "selected function help".


texts:
Using ArcGIS Pro 2.5.0 I've established a grid of points from a raster using the 'Raster to Point' tool and added the X and Y coordinates with 'Add Geometry Attributes'.
I now need to select a subset of points on a certain latitude. Using the 'Select Layer By Attribute' tool I specify the exact POINT_Y but the selection does not return results, in contrast to searching by OBJECTID. I assume this is due to rounding errors - the decimals displayed are not the decimals given when I click the actual value (in the case of the value in the image, I get -5.98330575 when I double click on the field).

I've tried both values and also tried to change the number of decimals displayed, but still selecting by POINT_Y fails. Is there another tool to achieve this or a way to set decimals when calculating the point coordinates in the first place?


texts:
I am working in a company where our team requirement is to do reverse geocoding lookups. Below is our requirement.

To give you more background, the 3 geolocation-related needs I have are as follows:

Lat/Long lookup: returning country, city, state, region, postal code
Address lookup: returning lat/long
Postal Code lookup: returning country, city, state, region, lat/long (I assume the lat/long of the postal code’s centroid is what is usually provided)

Which provider can I look into which will fulfill my above needs? I am looking for something which is used worldwide and are more reliable in terms of accuracy and coverage.


texts: Ok Button disabled until every input has been filled

texts:
I'm developping a QGIS Plugin where the User clicks on the map, my plugin make a buffer and propose a street name, a street side and a disctrict, then the user selects which administration he's from and the status of the building select. After this the User needs to puts 3 values to add to our database.

So, I have a box with several inputs from the user.

I would like to set the "ok" button enabled only when all the inputs are filled.

Is there anyway to 'check' if everything is filled in each function and tell the button "ok you can go". Like the Checkbox but without a check by the user?

Here under part of the code for the first line with an user's input (adpn):

        self.adpn_text = QLineEdit()
        self.adpn_text.textChanged.connect(self.get_values_adpn)
        self.layout.addRow(QLabel("Entrer l'ADPN (format = num):"), self.adpn_text)

    def get_values_adpn(self, text):
        self.adpn_values = text

    def click_ok (self, click):
        print("--- Récupère les données de l'User---")
        print('ADPN User', self.adpn_values)


texts: ArcGIS function to subtract cell from next cell in a raster?

texts:
Using ArcGIS 10.5.1 -
I have a raster with certain values (elevation in this case). I'd like to generate a raster that is based on each cell subtracted from its neighboring cell (the one to the direct left). Is there a function that would do this? Example:

Input raster:
 6    16    14
6    10    14
8    10    14

Desired output raster:

No Data   10    -2
No Data   4     4
No Data    2    4

Used Raster --> Float. Here is my header file


texts:
I'm having the same problem from this
previous question.

I keep getting some nulls after running expression. Not sure why it works on some dates

I have dates(strings) in the form 1/4/2005, 11/11/2004, etc

Here is code I used in field calculator:

btw data will be uploaded to django app using PostGeSQL

from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import QDate

@qgsfunction(args="auto", group='Python')
def parse_date_dmy(fromval, feature, parent):
    return QDate.fromString(fromval, 'dd/MM/yyyy')


texts: Using QGIS to create a table of correspondence between two area layers

texts:
I have two area layers over a map of Australia, and I want to know if it's possible to use QGIS to create a table of correspondences between them.

The Australian Bureau of Statistics provides tables of correspondences between areas, and lists the ratio of one area type that falls into another area type. The example below show the correspondence between Statistical Area 2 and federal electorates:

Is there a way to create a table like this using, say, the Field Calculator function in the Attribute Table once I've joined the two layers?


texts: Is it possible to open a raster or raster dataset contained in a file geodatabase with QGIS 3.X?

texts:
I know that it is possible to extract vector data from .gdb's. I have found some similar questions here and here or even here on StackExchange, but they were either not exactly what I was looking for or really old where the answer was, that QGIS cannot read raster data from file geodatabases.

It is 2019 now and 3.6.2 was just released. Is it now possible to open raster data from file geodatabases with QGIS? Or at least a way to check whether the geodatabase even contains raster data without having any ESRI products?


texts: How to upload my ArcGIS Server web map to free domain?

texts: i have created web map with arcgis server version 10.1, i just want to know how to upload this map from local url to free domain.


texts:
I'm trying to put a custom command into the Help menu on startup if it isn't already there. This is my code to do so:

private void OnCreate()
{
    // other stuff in OnCreate()...
    // registering my method for the right event
    m_appStatusEvents = m_application as IApplicationStatusEvents_Event;
    m_appStatusEvents.Initialized += new IApplicationStatusEvents_InitializedEventHandler(AddCmd);
}

private void AddCmd()
{
    UID myCmdId = new UIDClass();
    myCmdId.Value = myCmdIdValue;  // myCmdIdValue is my command's UID

    ICommandBars commandBars = m_application.Document.CommandBars;

    UID helpMenuId = new UIDClass();
    helpMenuId.Value = helpMenuIdValue; // helpMenuIdValue is esriArcMapUI.MxHelpMenu
    ICommandBar helpMenu = null;
    helpMenu = (ICommandBar)commandBars.Find(helpMenuId, false, false);
    // here, helpMenu.Count is 3
    if (helpMenu != null)
    {
        ICommandItem myCmd = helpMenu.Find(myCmdId, false);
        if (myCmd == null)
        {
            Object index = Type.Missing;
            helpMenu.Add(myCmdId, ref index);
            // here helpMenu.Count is 4, so my command was added
        }
    }
}

I think it's pretty straight forward: I'm looking for the Help menu and add my command if it isn't already there. In the debugger, I can see that an item is added to the Help menu, as its Count goes from 3 to 4.

Problem is: The command doesn't show up. I can add it manually, so I know it's loaded and correct, it just doesn't show up when I try to add it automatically.

Edit: I just tried to refresh the Help menu with (helpMenu as ICommandItem).Refresh() but that didn't help either.

Edit2: I just added this to debug:

for (int i = 0; i < helpMenu.Count; i++)
{
    Debug.WriteLine(helpMenu[i].Caption);
}

And it shows

&ArcGIS Desktop-Help
&ArcGIS Resource Center
&About ArcMap...
myCommand Caption

So it's definitely there, it just won't show!


texts: How can I create GeoJson map from OSM Data

texts:
I am trying to build a map of countries in the some countries in Europe (Spain/France/Portugal).

What is the best way to build such map:

Should I look for shp files, then concatenate them.
Is there any way to extract, and concatenate GeoJson of specific countries (with their boundaries data), then concatenate them.

It should be something similar to build such kind of maps. Any starting point ?


texts:
I'm trying to write a python script that will add a tif to the table of contents of my map.

(Eventually I want to convert to tool that

allows the user to type in the filename and then has the tool create the full path name, and
zooms to the extent of the newly-added tif, but baby-steps for now.)

My first stumbling block is a very basic one, seemingly.

I'm trying to add tifs, not layers, to my TOC. Do I need to use MakeFeatureLayer(), before I can use AddLayer() ?


texts:
I need to draw simple histograms in QGIS version 1.8 linked to data from an imported DB.

I need a double bar diagrams (eg: male/female - in/out) but I only find the possibility to set pie charts. Is it possible to draw histograms with high proportional to the value?


texts: Will ArcGIS Desktop support latest Python 2.x or Python 3.x versions?

texts:
I'm curious to know whether ArcGIS Desktop 10.0 will support newer versions of Python than 2.6. And for that matter newer versions of NumPy (ie, 1.4, 1.5 etc). I'm sure that this documentation is on the ESRI site somewhere, but I can't find it.

I have found that ArcGIS Desktop 10.0 will respect an existing Python 2.6 install as long as you have Numpy 1.3.0 installed (ref), and I'm guessing, at least until Numpy 1.5.0 is out and supporting Python 3, then no Python 3 for ArcGIS, but I'd still be interested to know if upgrading is possible when it's available.


texts: Exporting and clipping soil layer downloaded from ISRIC

texts:
How do I export and clip a soil layer downloaded from ISRIC?
I just downloaded a soil PH to ArcGIS, but I only need the 0-5cm mean layer, I found no place to export this layer to a raster/vector. And I only need the North American extent.
I am thinking 0-5 cm mean layer should be extracted first, then clipped to North America.


texts: Coding grass through spyder: "No module named 'grass.script'" error

texts:
I'm trying very hard to get spyder to run the grass module but have not succeeded.

I tried the grass documentation and these Q&As:

Using Spyder IDE to develop python for Grass GIS, get 'ImportError: No module named grass.script'
and
How do I connect my Python script external to the Grass GIS 7 program in Windows 10?

This is my code:

import os
import sys
import subprocess

grass7bin = r'C:\OSGeo4W64\bin\grass72.bat'
startcmd = [grass7bin, '--config', 'path']
try:
    p = subprocess.Popen(startcmd, shell=False,
                         stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out, err = p.communicate()
except OSError as error:
    sys.exit("ERROR: Cannot find GRASS GIS start script"
             " {cmd}: {error}".format(cmd=startcmd[0], error=error))
if p.returncode != 0:
    sys.exit("ERROR: Issues running GRASS GIS start script"
             " {cmd}: {error}"
             .format(cmd=' '.join(startcmd), error=err))
gisbase = str(out.decode('utf-8'))
gisbase = gisbase.rstrip()
os.environ['GISBASE'] = gisbase

#set up GRASS environment variables
grass_pydir = os.path.join(gisbase, "etc", "python")
sys.path.append(grass_pydir)

import grass.script.setup as gsetup

I ran spyder with python 3 and grass 7.2.0.


texts: How to correctly display trees in ArcScene?

texts:
The first image is a screenshot of the view in ArcScene.  I am trying to export this scene to a tif or jpg.  The scene is composed of the following layers:

A raster
Three polygon layers displayed as outlines only
One transparent polygon layer
A point layer displayed as trees based on height values within the point attribute table.

When I export the scene to a jpg or tif, the resulting error image is shown below with the gray overlay.  I was able to isolate the problem though--it appears the export works for everything except the points (displayed as trees).  If I export the scene with the points not displayed as tree, this also works.  Where am I going wrong?  How can I resolve the issue so that the trees are displayed correctly?


texts: Encoding issues (UTF-8 to latin1) in the production of shapefiles using QGIS API

texts:
We are facing encoding problems with shapefiles trying to convert UTF-8 encoding to latin1/ISO8859-1 using QGIS API.

In QGIS desktop over the layer added to the TOC doing "save as" and specifying the new encoding (latin1) produces a new shapefile with the encoding we expect, the dbf file is good since the text in attribute table is correct; although since we need a way to automate it over all the shapefiles, we have tried using API - quis.core and PyQt4.QtGui - for the processing, and it creates a new shapefile theoretically in latin1 encoding (properties says latin1 and .cpg file also), but when we inspected it, the text/characters are incorrect. We have used directly the QgsVectorFileWriter class, and seems it does not work properly over the dbf file..!?
Using QGIS 2.12 and python 2.7 (x86).

Any idea how to solve the weird behavior/issue?

Main code snippet:

app = QApplication([], True)
QgsApplication.setPrefixPath(r"C:\Program Files (x86)\QGIS Lyon\apps\qgis\\", True)
QgsApplication.initQgis()

layer = QgsVectorLayer(<path_to_input_shp>, "lyr_from_shp", "ogr")
if not layer.isValid():
    print ("Layer {} did not load".format(layer.name()))

writer = QgsVectorFileWriter.writeAsVectorFormat(layer, <path_to_output_shp>, "latin1", None, "ESRI Shapefile")
if writer == QgsVectorFileWriter.NoError:
    print ("success")

QgsApplication.exitQgis()


texts: MapTiler labels plotting too small in QGIS

texts:
I am using MapTiler vector tiles to show roads and suburb names in QGIS. If I export to pdf from a layout, the road names and suburb names show okay (albeit smaller than what is shown onscreen in QGIS). But if I plot directly to the plotter from QGIS (using the exact same QGIS layout), the road names and suburb names are so small that they can barely be seen, never mind being able to read them. Reducing the plot resolution in the QGIS layout increases the size of the road names and suburb names, but results in jagged text. I have experienced this issue in both 3.14 and 3.16.
I know that one possible solution is to display the basemap without the road names and suburb names and generate the names separately, but the fact that the road names and suburb names in the pdf file are significantly larger than those in the direct plot from QGIS suggests that QGIS is playing a role, not just the basemap itself.
The problem is similar to this one:
Label appears completely different in Layout-view than in exported image in QGIS using maptiler plugin


texts: messed up shapefile attribute table text

texts:
I have a .shp file from AutoCad imported to QGIS (using _MAPEXPORT > text in AutoCad) and this is what happens to text in attribute table - example:

{\Fsimplex|c0;O\Fsimplex|c238;Š\Fsimplex|c0; Dom\Fsimplex|c238;ž\Fsimplex|c0;ale
\H0.8x;Bistri\Fsimplex|c238;š\Fsimplex|c0;ka cesta 19
ID povezave 115006}

It seems it doesn't recognize that {\Fsimplex|... is a font definition. I went through some QGIS settings related to fonts and style but that doesn't seem to be a problem.


texts: Link multiple value from table to point shapefile

texts:
In QGis 2.18 is there a way to link (like join) a table with multiple value for each ID CODE with a Point shapefile with only 1 point for each ID CODE?? (picture in annex)

It would be fine also a temporary link in my machine memory to keep light the shapefile!


texts: Create new polygon from area bounded between two existing layers

texts:
I am trying to create a new polygon layer made up of the area of an existing polygon layer south of another layer that bisects it. This is a little difficult to explain, so here is a screenshot that I hope helps:

There is a river layer (gray) that bisects the county layer (blue). I want to create a new polygon in QGIS that includes all of the county area south of the river layer. Both layers are polygons.


texts: "Geometry type does not exist" on a Postgis installation on Fedora 18

texts:
Just installed PostgreSql and PostGIS over a brand new installation of Fedora 18. I did it this way:

sudo yum install postgresql postgresql-server postgrsql-contrib postgresql-devel pgadmin3

Once installed PostgreSql, I created the postgres user as usual.

Then, I installed PostGIS:

sudo yum install proj-devel gdal postgis

As a last step, I've created a new user (mydbuser) and a new db owned by that user:

psql -U postgres -W
CREATE DATABASE mydb OWNER mydbuser;
\q
createlang plpgsql mydb

Now, I've tried to execute the following line to use PostGIS on mydb:

psql -d mydb -U mydbuser -W -f /usr/share/pgsql/contrib/postgis-64.sql

but I get a long list of errors of the same type, this is the first line:

psql:/usr/share/pgsql/contrib/postgis-64.sql:7751: ERRORE:  il tipo "geometry" non esiste

In English, the message is 'ERROR: "geometry" type does not exist'.
I've also tried to execute the code as the postgres user, but nothing changed. Any ideas?

Diego


texts: Optimizing projection for accuracy in calculating distances

texts:
I receive a set of LatLong points and I want to optimize the projection to Euclidean points so as to try to keep (distance in Euclidean space)/(greater circle distance in latLong) with variance as low as possible.

The Mercator projection error in distances (east/west stretching) increases as you get farther from the equator. Is it common to calculate the projection with a different latitude than then equator as the "center"? If yes any methods/material on this?

The set of LatLong points has varying latitudes but always has a |minLat - maxLat| < 10

Usually the radius of a circle containing the points is 100km and not exceeding 10,000km. I do not expect the projection/process to work as well when the points are more spread out.


texts:
I'm trying to use locate function from Leaflet in R, I'd like a marker to show up after clicking on locate button but the following code only centers and zooms the view on the user's geolocation. I'm new to JS, and I thought passing the drawMarkers: true, from the manual  would map markers:
library(sf)
library(dplyr)
library(leaflet)
library(htmlwidgets)
library(htmltools)
library(leaflet.extras)

punto = st_sfc(st_point(c(-96.951239, 19.4600523))) %>% st_set_crs(4326)
punto_b = punto %>% st_transform(32614) %>% st_buffer(50000) %>% st_transform(4326)  %>% st_as_sf

leafIcons <- icons(
  iconUrl = "https://gitlab.com/datamarindo/anps/-/raw/master/logo_prona.png",
  iconWidth = 100, iconHeight = 100)

l <- leaflet(data = punto) %>% addTiles() %>%
      addMarkers(dat = punto, icon = leafIcons) %>%
      addPolygons(data = punto_b, fillColor = "yellow", weight = 4) %>%
      addMeasure() %>%
      addMiniMap() %>%
      addEasyButton(easyButton(
         icon="fa-crosshairs", title="Locate Me",
         onClick=JS("function(btn, map){ map.locate({setView: true, drawMarker: true});}")))


texts: In QGIS Web Client - How to click on a photo in an attribute table and load into a new browser tab?

texts:
I have a map loaded into my QGIS web client, and have images loaded into my attribute table for a specific layer. Works great - when I click info on one of the layer attributes, it shows the image.

BUT - the image is very small. I would like to be able to click the smaller version of the image and have it load the original image into a new browser tab.

I am able to right click and click "open image in new tab" and it works, but I'd like anyone who looks at the map to be able to just click it and automatically have the image pop up in a new window.

Is this possible?


texts: Removing values from shapefile with ArcPy?

texts:
I have a large shapefile with lots of duplicate records. However, each duplicate record has an accuracy rating. I want to remove all duplicate records with a lower accuracy rating.

This would be simple but there can be duplicate records with the same accuracy rating and I want to leave them in. This might be easier to show than to explain:

Key Accuracy

ABC 100

ABC 100

ABC 50

XYZ 100

XYZ 75

XYZ 75

123 50

123 20

123 20

In the above table I would want to keep the top two ABC records as they both have the best accuracy, but remove the last one. I would want to keep the top XYZ record but remove the other two. And I would want to keep the top 123 record as, although it doesn't have the best possible accuracy, it has the highest accuracy for that batch of duplicate records.

How do I do this using ArcPy?


texts: Adding a point shapefile to a raster image with no spatial information

texts:
I have a raster image (JPG) of a fictional map from a fantasy book.

I want to add a point shapefile in order to show main events from the book upon the map.

My map has no spatial information (as it is fictional) and therefore I am unable to "start editing" the image, which prevents me from adding points.

Any advice to how I can add a points to a JPG image with no spatial information?

This is the map I'm trying to add a point shapefile to.

http://www.ursulakleguin.com/EarthseaMaps/


texts: pyQGIS 2.0: Sending signal when identify button is pressed via gui

texts:
I have a question regarding the qgis identify tool for version 2.0. I would like to know if it's possible to send a signal when the actionIdentify is triggered by the user via the gui.

I'm trying to call a python function within a plugin if this action is triggered.

I've tried two different ways so far with no success:

1)

def initGui(self):
    ...
    self.iface.actionIdentify().triggered()
    self.idTool = self.canvas.mapTool()

def myFunction(self):
    if self.canvas.mapTool() == self.idTool:
        QgsMessageBox.information(None, 'Info', 'Tool is IdentifyTool')
        self.anotherFunction()

def anotherFunction(self):
    QgsMessageBox.information(None, 'Info', 'Successful connection')

This gives no error, but the QgsMessageBoxes don't show up, thus I am not successfully connecting to the function. It's also problematic, because my way of setting the idTool parameter is by first triggering the tool. This means that even if everything runs as I thought it would, the if statement would always be true, as the self.canvas.mapTool() would surely be the idTool.

this gives me a PyQt4.QtGui.QAction object:

idTool = qgis.utils.iface.actionIdentify()

this gives me a qgis.gui.QgsMapToolIdentify object:

qgis.utils.iface.actionIdentify().trigger()
anotherIDTool = qgis.utils.iface.mapCanvas().mapTool()

I can't successfully compare the former with my current map tool as returned from mapCanvas().mapTool(). However, I don't want to have to trigger the tool first to get the proper object.

2)

def __init__(self, iface):
    ...
    self.idTool = QgsMapToolIdentify(self.canvas)
    ...

def initGui(self):
    ...
    self.idTool.mapToolSet.connect(self.myFunction)

def myFunction(self):
    currentTool = self.canvas.mapTool()
    if currentTool == self.idTool:
        self.anotherFunction

My problem here is that the two QgsMapToolIdentify parameters point to different object numbers, so I never get past my if statement. What am I missing about the way these QgsMapToolIdentify objects are named?

UPDATE I am revisiting this issue for QGIS 2.8. Does anyone have any ideas?


texts: Why are there no line breaks in my labels even though I added \n?

texts: I'm having an issue with creating multi-line labels in QGIS.  I believe I have the syntax right in the label expression because the output preview looks correct.  However, on the map, everything displays on the same line.  If anyone can help, let me know - I'm new to both QGIS and StackExchange.  Screenshots are below.


texts:
Is it possible to get the decorator scale bar onto the map canvas inside a standalone PyQGIS application? It is possible to get a scale bar on the QGIS Desktop version by View -> Decorations -> Scale Bar
I saw this post: Access / modify QGIS copyright decoration through python but this is solely only for the copyright decorator. Was wondering if there is anything similar for the scale bar.


texts: How to check if shapefile is empty using ArcPy?

texts:
I have a number of shapefiles in a folder but some of them are empty i.e. there are no rows within the shapefile.

Is there a simple way to find those empty shapefiles using ArcPy?


texts: Add polygon (buffered line) created by turf.js to map in Leaflet?

texts:
I am trying to create a buffered line layer by creating GeoJSON line string from  routing line and buffering it with turf.js and then show it on the map.
Here is my code:
route.on('routeselected', function(e) {
    var r = e.route;
        var lineCoordinates = [],
            i,
            latLng;
        for (i = 0; i < r.coordinates.length; i++) {
            latLng = L.latLng(r.coordinates[i]);
            lineCoordinates.push([latLng.lng, latLng.lat]);
        }

        let lineString = L.polyline(lineCoordinates,{color: 'black'});
        var buffered = turf.buffer(lineString, 500, {units: 'miles'});
        var bufferedLayer = L.geoJSON(null);
        bufferedLayer.addData(buffered);
        bufferedLayer.addTo(map);
        });

Any suggestion?


texts: QGIS 'Save As' to Spatialite is changing integers into strings

texts: Using QGIS 2.14.0, Windows 10, when I use 'Save As' to save a polygon shapefile as a Spatialite layer (all options left as default), all of the input integer fields are incorrectly output as text string fields.  This sure seems like an error to me.  I reviewed the Spatialite web site, but didn't see anything.  Am I doing something wrong?


texts: Difference between maxBufferCount and maxRecordCount ArcGIS

texts:
I have defined featureLayers for my web map application, However it is very slow loading all polygons of mapLayers.

would increasing  "maxBufferCount" or "maxRecordCount" be a partial solution?

What is the difference between these two?


texts: QGIS make lines with offsets of y and azimuth from source and target coordinates

texts:
I'm very new with QGIS and currently stuck with such problem in my work. Maybe some one knows how to help me?

I have made two separate layers.

Database layer with sites and antennas, containing site center points and azimuths with antennas. So I have created symbols representing such antennas by creating an offset from y points by the size of the symbol in a direction of azimuth data. It looks like this:

I need to analyse traffic between these antennas, which has enormous amount of data, so in order to not overload the system with data I import some data of wanted object (a new layer in QGIS) containing traffic destination with latitudes of source antenna and target antenna coordinates. together with that I import the azimuth info of both points.

So I create the lines using the function:

It looks like that when I enable both layers on the map:

So where is the problem?

If you look at this lines I cannot see from which sectors to which sectors of another object the line goes to because they all are drown from the centres. How can I make some offsets in the direction of known azimuth in the source point and target point?

In another words:

have lines like this: point1xy....point2xy

I want:

point1xy+y offset (-10) in direction (azimuth1).....point2xy+y offset (-10) in direction (azimuth2)


texts:
I'm trying to make a spatial query, but I think something is wrong with it:

I have two tables, places (point) and cities (polygon)

When I try to get the points within 500m to one point, there is no problem with st_dwithin:

select * from places
where st_dwithin(geo, (select geo from places where id=2), 500)

But what If I want to get the polygons within 500m of one polygon, and then, all the points inside this selection?

I'm trying to do the first part: Polygons within 500m of my given polygon:

select * from munisp
where st_dwithin(geo, (select geo from munisp where id=275), 100)

And I get the error:

  BOOM! Could not generate outside point!

The 'geo' column are Geography data types, in both tables, Polygon, 4326 in one, and Point, 4326 in the other.

Am I missing something?
Any idea about how could I get the points inside this, and neighbour polygons?

The point table was geography since the beginning, but the polygons one was cast from geometry (created geography column 4326, and filled with update polytable set geo=the_geom::geography)

Could that be the problem?

Having tried the query:

SELECT p.* FROM munisp m, munisp p WHERE m.id=275 AND ST_DWithin(m.geo, p.geo, 100);

I get the same error.

  BOOM! Could not generate outside point!

And when I try:

SELECT p.* FROM munisp m, munisp p WHERE m.id=275 AND st_touches(m.geom, p.geom);

*Note the 'geom' column, not 'geo'. I get 0 rows. What means that no polygon is touching #275, so what is false, and it should return multiple polygons.

Could it be due to an error in the data?


texts:
How to receive GPRS data and view it on a webpage?

I am doing a GPRS Datalogger. It sends Data of some sensors to a web server using GPRS. I want to be able to view this data on a webpage in a browser say Google Chrome. Also I need the data which is received by server to be stored in a database. How should I proceed with the designing of Front end for my application? Should I use HTML 5, PHP and MySQL? Is there an example code which receives data from GPRS and shows it on a webpage?


texts:
I have a Linestring layer for which I have build a QgsSpatialIndex with the following code:

my_linestrings_lyr_idx = QgsSpatialindex()
for f in my_linetraings_lyr.getFeatures():
    my_linestrings_lyr_idx.insertFeature(f)

Now I want to find the nearest neighbor for a QgsPoint that is definitly located on one of the linestrings layers features, code:

#assume pt is of type QgsPoint
fid = my_linestrings_lyr_index.nearestNeighbor(pt, 1)
r = QgsFeatureRequest(fid[0])
nearest_feat = my_linestrings_lyr.getFeatures(r).next()

I expect nearest_feat to be the feature on which pt lies in any case (because this is definitly the nearest neighbor to pt), but in some cases that I can't exactly specify, nearest_feat does contain other features, which does not contain pt.

Do I use the nearestNeighbor() method incorrectly in some way?


texts:
I'm having trouble locating any discussion or comparative benchmarking of different raster file formats (e.g., for use in data analysis in R). Does anybody have any insight into why particular formats may be faster or slower? Or should the differences be minimal?

Specifically, I'm interested in if converting a raster (e.g., a GEOTIFF file) to a different format (e.g, a netCDF) is ever worthwhile for the purpose of speeding read/write and other operations.


texts: How to count zero values in a table in rows and write it to new field?

texts:
This is a very simple task to do but I can't understand the right syntax.

I have a shapefile whose attributes are similar to the following:

FID Shape   FIELD1 FIELD2 FIELD3 ...
0   Polygon 0      1      0
1   Polygon 3      0      7
2   Polygon 3      4      7
...

The number of fields and their names are always different.

I need to create a new field (let's name it NUM), and count the number of zeros in each row.

Example output:

FID Shape   FIELD1 FIELD2 FIELD3 NUM
0   Polygon 0      1      0      2
1   Polygon 3      0      7      1
2   Polygon 3      4      7      0

I know how to create a new field, however I am not clear on the next steps.

The working code:

#path is path to shape file
def a(path):
  fields = arcpy.ListFields(path,"FID_*") #FID_* is wildcard to select a fields name
  arcpy.AddField_management(path, "NUM", "SHORT") #create a field with name NUM
  cursor= arcpy.UpdateCursor(path)
  for row in cursor:
    count=0
    for field in fields:
      a= row.getValue(field.name) #take a value
      if a==0: #if value=0 then value=value+1
        count+=1
    row.setValue("NUM", count)
    cursor.updateRow(row)
  del row
  del cursor

Thanks blah238, now I can eat pythons!


texts: Canceling task currently running in background using PyQGis?

texts:
I am trying to close QgsProject on loading layer (if password and username from QgSProject - datasource are  incorrect) before QgsCredentialDialog is shown.

on signal :

  QgsProject.instance().loadingLayer.connect(self.load_log)

If password and user are correct -> pass
if password and user are not  correct: close the project

I am using :

    project = QgsProject.instance()
    project.clear()

but QgsCredentialDialog still shows up.

after clear project if I use :

  iface.actionExit().trigger()

This shows up:

If I click yes, QgsCredentialDialog is shown again

How do I close the project and everything connected with that project (signals,actions,tasks,connections,...)?

NOTE:

QgsTaskManager().cancelAll() # NOT WORKING
tasks = QgsApplication.taskManager().activeTasks()
for task in tasks:
    print (task.canCancel()) # False
    print (task.description()) # Loading “C:/Users/User/Desktop/some_project.qgs”
    print (task.flags()) # <qgis._core.QgsTask.Flags object at 0x000002330527E4C8>

    task.cancel() # NOT WORKING

Since task.canCancel() is FALSE (I think that is main problem here), is it possible to change task.canCancel() into TRUE?


texts:
I have an ArcMap 10.3.1 basic License and therefore cannot export data to SDE's connections.  Is there an open source tool that I can use to export shapefiles to an sde?  I have found shp2sde and sde2shp but these tool are outdated because of this statement in the documentation:

  Note: This function is not supported on 64-bit Windows, Linux, or HP-Itanium.

I'm using Windows.

If there is not a tool do I have any alternative other than upgrading to Advanced?

Is there anyway to do this without upgrading my existing license.  i.e some light weight opensource tool


texts: How to permanently change attribute values using PyQGIS?

texts:
I'm trying to iterate through all the features of a layer to change the value of some of its attributes. I'm using the following code:

from qgis.core import *
import os

working_folder = r'C:\myFolder'
fc_Opt_Prospect = os.path.join(working_folder, 'Equipm.gdb|layername=Prosp')

layer = QgsVectorLayer(fc_Opt_Prospect, 'Optimity_Prospect', 'ogr')

layer.startEditing()

for f in layer.getFeatures():
    layer.changeAttributeValue(f.id(), 4, 3)
    layer.changeAttributeValue(f.id(), 5, 6)

layer.commitChanges()
QgsProject.instance().addMapLayer(layer)

When I check the attribute table of the layer I can see that no value has been changed. No error is thrown by the script.


texts:
How can I Create a spatial binning similar to or directly using Uber H3 from QGIS UI?
Are there plugins available to achieve the desired task?


texts: Merging Raster data with displaced points vector data

texts: I have two data-sets.  One is a point vector of which for respondent confidentiality reasons the points are displaced by 5KM. And then I have the satellite raster data. To account for the displacement error in the vector data I wanted to aggregate the raster data to 5KM resolution before merging the two data-sets. How do I aggregate the raster to 5KM resolution in QGIS?


texts: Creating point shapefile from Excel with Fiona

texts:
I am attempting to use fiona and xlrd to convert a complicated excel spreadsheet to a shapefile (point).  I am trying to figure out how, during the .write call, to send the list of coordinates I've pulled out of the spreadsheet and use them to draw the geometry.  Here's an abbreviation of what I've got so far:

shapefile_dir = r'C:\Temp'
point_shapefile = os.path.join(shapefile_dir, '{}_export_p.shp'.format(datetime.now().strftime('%Y-%m-%d %H:%m')))

schema_properties = {
    'guid': 'str:50'
    'dateTime': 'str:16'
}

point_schema = {'geometry': 'Point', 'properties': schema_properties}

with fiona.open(point_shapefile, 'w', crs=from_epsg(4326), driver='ESRI Shapefile', schema=point_schema) as layer:
    with xlrd.open_workbook(wkbk) as wb:
        cs = wb.sheet_by_name('Master')
        keys = [cs.cell(0, col_index).value for col_index in xrange(cs.ncols)]
        dict_list = []
        for row_index in xrange(1, cs.nrows):
            d = {keys[col_index]: cs.cell(row_index, col_index).value for col_index in xrange(cs.ncols)}
            dict_list.append(d)

        for idict in dict_list:
            try:
                guid = idict['A']
                dateTime = idict['DATE']
                declat = idict['DECLAT']
                declon = idict['DECLON']
                try:
                    layer.write({'geometry': (declat, declon), 'properties': {
                        'guid': guid,
                        'dateTime': dateTime
                except: #except clause for writing the geometry and attributes
            except: #except clause for getting values out of dictionary list

I'm a little overboard with xlrd because I need this to scale to read a workbook with n number of sheets.

But currently I'm getting an error

tuple indices must be integers, not str

Which of course means I'm missing something very basic.

My values are being fetched from the excel spreadsheet correctly, so I think I'm almost there.  But I'm not sure I'm thinking through this procedure correctly.  Should I be stashing the xlrd operations inside of the fiona write loop?  How do I get from that list of XYs, cast them to float and then write them as a new row in the shapefile?

I reviewed Fiona - Preffered method for defining a schema

but I guess I'm not clear on how to adapt some of the methods in there to my own data.


texts: I would like to insert lines into a table that I have in my PostgreSQL db from various tables in my computer

texts:
I have a table on my PostgreSQL database that I would like to update by inserting the missing lines that I have in many other shapefile tables. I was able to connect to my server PostgreSQL from the Windows' shell (psql -h localhost  -U quest  -d "test" -p 5432) and import some shapefiles to update my table (to insert the data from pgAdmin), but instead of importing all my tables into PostgreSQL, I would like to update my table directly from the shell. The problem is that I do not know how to go about it. Specifically, I don't know how to tell PostgreSQL the source of my tables in the "from".
In pgAdmin I use :
INSERT INTO bal.bal (bal_pk_uid, bal_agg_id, geom)
            SELECT  bal_pk_uid, bal_agg_id, geom
    FROM bal.bal_2 bal2
    WHERE bal2.bal_pk_uid is not null and not exists (select 1 from bal.bal balglobale where balglobale.bal_pk_uid = bal2.bal_pk_uid);

but in the shell how am I supposed to tell the path to my source "C:"CARTO"\bal.shp bal2"?


texts: Calculating share of overlapping areas via aggregate function in a virtual field

texts:
I'm trying to calculate the shared area for several overlapping areas. I need a dynamic solution in order to get automatically calculated statistics when drawing new polygons.

The target area is in this case the layer which contains the zoning plans. The other layer specifies the type of landuse (Markanvändning). The ultimate goal is to have a virtual field in the zoning layer, which grabs the information from the landuse-layer to define the share of different landuse-types (area) within the plan.

((aggregate(layer:='Markanvändning_e90475ff_88be_43f2_8eae_b14419c74d08', aggregate:='sum', expression:= (area(intersection($geometry, geometry(@parent)))), concatenator:=',',  filter:= "kategori" IS 'Offentlig'))/ "area_DP")*100

My approach was to create a virtual field within the zoning plan layer with an aggregate function, which calculates the share with the help of an intersect function. I created a test-layer and everything works fine if I just have one zoning-plan, but when I draw several plans I just get a result for the first row in my zoning plan layer.

What is going wrong here?


texts: I need to use traverse method in my program on visual studio(WPF) so that the user can get their target point coordinate by using their initial data. In my knowledge there is two type of traverse that is open traverse and closed traverse. I search on the internet there mostly about the closed traverse. Although the open traverse cant be check for error i still need to use it for my project requirement. Can someone share me a simple open traverse calculation that i can use for my reference and program?


texts: New to QGIS here.  I'm trying to figure out how to add information to set points.  I'm trying to use them as leads, so I can put lead name, address, phone number, whether they've been talked to, etc. How exactly would I go about doing that?


texts:
.osm.pbf it`s the extension for OpenStreetMap extract that I get from a site that takes daily extracts from osm: http://download.geofabrik.de/

How to import these files in QGIS? On Linux or Windows? The more simplest variant, if it exists.


texts: Is there a "clone view" functionality in Qgis like in MapInfo

texts:
I am currently working in Qgis (2.18.11) and I would like to view mi map window duplicated so I can represent a thematic map with different values in each window and compare both.
I usually do this in MapInfo (Map > Clone view). See the example:


texts: QGIS: How to convert date values with format d-m-yyyy in a string field, into date values with format yyyy-mm-dd in a date field?

texts:
I have date values with the format d-m-yyyy in QGIS.
They are in string field.
I want to convert them into date values with the format yyyy-mm-dd.
I want them in a date field.

Can this be done without doing something in Python?