texts: Error when looping Polygon to Raster with python through different fields using ArcGIS Desktop

texts:
I've been trying to create a python script to loop through numerous shapefiles converting various fields to raster outputs.

I've been working off a lot of what has been said in this previous thread on the same issue -> Python loop Polygon to raster not working?

But can still not get it to work. I keep getting the dreaded 999999 error.

The code I've been trying is shown below:

#import arcpy module
import arcpy
from arcpy import env

#set working environment
env.workspace = "C:\\Users\\s175850\\Documents\\GIS\\GIS_Data\\Source_Methodology\\Pesticides\\Application_Term\\2008"
Dir = env.workspace

#List feature classes
fclist = arcpy.ListFeatureClasses()

#loop
for fc in fclist:

  output = Dir + "\\r" + fc + "_ras"
  cell = 5

  #process polygon to raster
  arcpy.PolygonToRaster_conversion(fc, "January", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "February", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "March", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "April", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "May", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "June", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "July", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "August", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "September", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "October", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "November", output, "CELL_CENTER", "", cell)
  arcpy.PolygonToRaster_conversion(fc, "December", output, "CELL_CENTER", "", cell)

print "finished polygon to raster"

I'm not sure if there would be issues with outlining the 12 polygon to raster conversions for the 12 fields like I have done so I tried it with just one polygon to raster run for one field name. I still managed to get the same error.

I am using ArcGIS 10.


texts: R gstat krige() NA predictions increasing with larger nmax

texts:
I've been playing around with a King County house dataset posted on Kaggle (link: https://www.kaggle.com/harlfoxem/housesalesprediction) and tried to predict house prices with Kriging. The entire coding is shown below:

library(gstat)
library(sp)
library(ggplot2)
setwd("~/Documentos/R/Kriging")
dados = read.csv('kc_house_data.csv')

dados$price_sqft = dados$price / dados$sqft_living

ggplot(dados, aes(x = long, y = lat, size = price_sqft)) +
  geom_point(alpha = 0.2) + scale_size_continuous(range = c(0.5, 5))

set.seed(1337)

treino_ind <- sample(seq_len(nrow(dados)), size = 0.80*nrow(dados))
dados_treino = dados[treino_ind,]
dados_teste = dados[-treino_ind,]

coordinates(dados_treino)  = ~ long + lat
coordinates(dados_teste) = ~ long + lat

vrg = variogram(log(price_sqft)~1, (dados_treino))

vrg

vrg.fit = fit.variogram(vrg, model = vgm(1, "Sph", 1, 1))
vrg.fit

plot(vrg, vrg.fit)

vrg.krige = krige(log(price_sqft)~1, dados_treino, dados_teste,
      model = vrg.fit, nmax = 1)

At first, I ran the code without the nmax = 1 bit, but saw some posts where the inclusion of some restriction on the amount of calculations was advised and finally added it. The thing with this is that if I increase the nmax parameter, kringe() function returns many NA. For instance, when nmax = 3, there was 216 NAs. With nmax = 10, there was 766 NAs. With nmax = 50, there was 2568 NAs (to give a little perspective, the entire test set is composed by about 4300 obs).

My question is: does that actually make sense or is it a bug?


texts: FeatureClassToFeatureClass_conversion where_clause ignored in PyCharm console

texts:
When I call arcpy.FeatureClassToFeatureClass_conversion from the PyCharm python console, the where_clause argument is ignored. When a value is supplied, it does not return an error message, nor does it apply the where_clause to the operation. For example:
In[4]: arcpy.FeatureClassToFeatureClass_conversion(r'C:\Workspace\test.gdb\wells', r'C:\Workspace\test.gdb', 'test1', where_clause='This is not sql')
Out[4]: <Result 'C:\\Workspace\\test.gdb\\test1'>
Stranger still, this command works correctly (in the above case, it returns an error) when executed from the Windows command interface (but using the same python environment). The where clause is only ignored in the PyCharm console. I am using arcpy from ArcGIS Pro 2.5.0, and PyCharm 2020.2.
I don't understand why this is happening, and it is compromising some of my workflows.
Can anyone replicate this error, or suggest a solution?


texts: Add WMS in QGIS gives an error with Map request failed [error:Host http not found url:http://http//

texts:
I got an URL from a friend who is working with Oracle MapBuilder. http://[IP]:11188/mapviewer/wms
I tried to get the service via ArcMap and it was working.

Then I tried my beloved QGIS and there is a problem. I get the connection running and all layers are shown in the "Add Layer(s) from a WM(T)S Server". I choose a layer and it gets integrated in the Layers window. But somehow there is no map. It is the same problem for all layers.

When I check the log Messages I can see the problem

Map request failed [error:Host http not found
url:http://http//[IP]/:11188/mapviewer/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&BBOX=35.02831983698249729,31.49752950471350132,35.04319025529049725,31.50711888517620096&SRS=EPSG:4326&WIDTH=982&HEIGHT=633&LAYERS=GT_BLOCKS_MB&STYLES=&FORMAT=image/png&DPI=120&MAP_RESOLUTION=120&FORMAT_OPTIONS=dpi:120&TRANSPARENT=TRUE]

The beginning of the URL is the problem

http://http//

and the part with the port

/:11188

.

How can I fix the problem???


texts: Field Calculation Python acceleration of gravity?

texts:
In my problem i have a column DISTANCE and a column TIME and i want to calculate with python in field calculator how much time it needs to make that distance (for example 1000m) if you are moving with the acceleration of gravity.

Since the field calculator does not natively support the square root function, I have no idea what to do.


texts: I have two shapefiles, one polygon one points.  The points represent specific field sampling locations.  The polygons represent land divisions.  The land divisions have a unique 16 digit identifier.  The plots have numbers unique to each polygon.  Ie:  polygon 0401010011190084 has plots 1-4.  I need to create a field in the plots shapefile that will be populated with the polygon number in which it lies.  Typing them in by hand sucks.  I need to populate this field in order to get the plots to load correctly into a database.  How can do this easily?    The plots shapefile needs to remain as points only.


texts: I have a polygon layer with multiple sectors, I created for each one a pie chart. My problem is some of thoses pies are not displayed, but when I zoom I find them. I would be pleased if you tell me how to have them all displayed.


texts: i can't find a way to export quadkey tiles using Qgis 3 and 2. I can export it to XYZ tiles and MBtiles, but i need quadkey tiles for MSFS2020


texts: Icon missing when adding custom combobox command to Toolbar

texts:
I have a custom command:

class ATL_NO_VTABLE CSliceSelector :
    public CComObjectRootEx<CComSingleThreadModel>,
    public CComCoClass<CSliceSelector, &CLSID_SliceSelector>,
    public ISliceSelector,
    public ICommand,
    public IComboBox
{
   ...
}

I also have a custom Toolbar where I try to add my custom command:

class ATL_NO_VTABLE CSeismicToolbar :
    public CComObjectRootEx<CComSingleThreadModel>,
    public CComCoClass<CSeismicToolbar, &CLSID_SeismicToolbar>,
    public ISeismicToolbar,
    public IToolBarDef
{
   ...
   STDMETHOD(GetItemInfo)(long pos, IItemDef * itemDef) {
        if (itemDef == NULL)
            return E_POINTER;

        if (pos == 0) { // CSliceSelector
            itemDef->put_Group (VARIANT_FALSE);
            HRESULT hr = itemDef->put_ID(CComBSTR(L"{96A2B223-BF0D-410C-B22B-892FED3D9A15}"));
            if (FAILED (hr)) {
                return hr;
            }
        }
      }
}

However when I start ArcScene; I see a red circle icon and the text "Missing": , instead of my command :-(.

However if I add it manually from Customize toolbars; it works fine.

I have tried both returning a null pointer and a valid bitmap:

CSliceSelector()
{
     m_hBitmap = ::LoadBitmap(_AtlBaseModule.GetModuleInstance(),MAKEINTRESOURCE(IDB_SLICE));
}
STDMETHODIMP CSliceSelector::get_Bitmap(OLE_HANDLE * Bitmap)
{
     if (Bitmap == NULL)
        return E_POINTER;
   *Bitmap = (OLE_HANDLE) m_hBitmap;
     return S_OK;
}

I have other commands that are not combobox; and they work fine.


texts: Adding file name of mapinfo file into its attribute table using MapBasic?

texts:
I have a task where I need to take the file name (no extension) of a MapInfo file and place it into a new field within its attribute table.

I would normally use Arcpy to do something similar, but I am unfamiliar with MapInfo and have a basic understanding of QGIS.

The reason behind this task is because we will merge all the MapInfo files and will need a reference to each of the features. We have a few hundred of these files so automation is preferred.


texts: Displaying raster Background Value as hollow for more than one set of RGB values using ArcMap?

texts:
I have two overlapping .ecw layers of the same area, though one is clipped in the middle, with a white edge.
When I set the background value to 255,255,255, and define it as hollow, it only hollows that specific RGB value, when in reality, there are other "white" values in my layer, such as:

254,254,254; 253,253,253.. all the way to 249,249,249.

Is there a way to define the display background value to include more that one unique RGB set of values?

I know I can work around this (clipping, raster calculator, etc.), but I want to retain only one original .ecw file without producing duplicates, therefore I only want to use Arcmap symbology or image rendering tools.

What I need is exactly like the Custom transparency options in QGIS, but for ArcMap.


texts: How to get the bounding box (MultiPolygon) of a geographic region?

texts:
I need to execute a BigQuery database query whereby I'm trying to ascertain if a given lat, lon coordinate is within one of a few US counties.

One of those counties is Collin County in Texas, shown here in Google Maps:

BigQuery supports GIS functions I'd need like ST_CONTAINS, but I can't figure out how to get, say, a GeoJSON representation of that polygon Google Maps is showing from either that service or OpenStreetMap.

The closest I've gotten with OpenStreetMap + Overpass is the following query:

[out:json][timeout:180];
relation(1837688);
out geom;

But the polygon I try to create from that response using members doesn't line up with what OpenStreetView has https://www.openstreetmap.org/relation/1837688.

Is there a more general solution I can use to create these bounding polygons for a region?


texts: Getting Negative Topographic Wetness Index (TWI) values in SAGA GIS?

texts:
I am getting negative topographic index wetness values in SAGA GIS - anyone else every get this?

Not all negative values, but the minimum of the TWI range is negative for some steep watersheds we are examining.  Ranges get as negative as (approx.) -1.5.

For this to happen, the value resulting from (specific catchment area / tan slope ) must be less than 1, since the ln of values between 0 than 1 = a negative value.

It is a very steep watershed, so we have slopes that approach 1.57 radians (almost 90 degrees).  In fact, the maximum slope for the watershed in question is 1.509 radians (or 86.459 degrees).  And the minimum Specific catchment Area value is 2.2175m.  So we can see that if we calculate ln(2.2175/tan 1.509), which further simplified is ln(2.2175/16.16159) and then ln(0.137208), we get a negative value (-1.8 or something) that shows the theoretical absolute minimum TI value for our watershed.  Its theoretical cause the steepest slope and smallest SCA don't combine, but some set of close values do, to get a real minimum of approx. -1.5.

To get more specific, in our steep watershed, the DEM has a resolution of 8.86 m and the Total Catchment Area minimum calculated by SAGA is this cell size squared (78.4996 sq m).  Which makes sense, as there are cells in the catchment that don't accept flow from any other cells and therefore have a Total Catchment Area equal to the cell area.  Now, for these cells, at localized peaks in the watershed (or even maximum peaks at the edges), the Specific Catchment Area should be equal to the Total Catchment Area for these cells divided by the contour width of downslope flow paths.  These peaks flow to all neighboring cells, this 78.4996 is divided by (as currently coded in SAGA) 35.4 (which is 8 * 0.5 * 8.86).  This gives the Specific Catchment Area minimum of 2.2175m.

My main problem with this negative TI result, is that I don't think it is possible, or at least have never really seen negative TI values talked about/graphed in a paper.

Any thoughts form anyone?


texts:
I have a point and a polyline. I need to find the perpendicular distance between each segment of the polyline and the point. By each segment of the polyline I mean that, each polyline may contain more than one line and I need to find every line of the polyline and find the perpendicular distance. I know its really easy using ICurve.QueryPointAndDistance Method. But to use this method, I need to create polylines from each segment of the actual polyline. But it is really time consuming for more than 13k polyline features.

So I tried to create line from every segment of the polyline. And calculate perpendicular distance from a point to the line. but this distance is not same as the distance I get from ICurve.QueryPointAndDistance Method's distanceFromCurve.

Question:
Is there any better option to calculate perpendicular distance that will be same as the ICurve.QueryPointAndDistance Method's distanceFromCurve? Or what is the reason I am not getting the same result as ICurve.QueryPointAndDistance Method's distanceFromCurve?

I am using arcObjects and c#. the input feature layer is always projected. So I am getting X, Y value in meter/feet. For calculating perpendicular distance I am using following method:

private double GetPerpendicularDistance(Line line, Point point)
{
    var a = line.Point1.Y - line.Point2.Y;
    var b = line.Point2.X - line.Point1.X;
    var c = (line.Point1.X * line.Point2.Y) - (line.Point2.X * line.Point1.Y);

    var x = point.X;
    var y = point.Y;

    var d1 = Math.Abs((a * x) + (b * y) + c);
    var d2 = Math.Sqrt((a * a) + (b * b));
    var d = d1 / d2;
    return d;
}

Line is an internal class.

internal class Line
{
    internal Point Point1 { set; get; }
    internal Point Point2 { set; get; }
}

Point is another internal class.

internal class Point
{
    internal double X { set; get; }
    internal double Y { set; get; }
}


texts: Loading *.adf files into QGIS using Python

texts: I am trying to write a loop to load .adf files into QGIS and create a contour map from each one. The response to an older question about loading .adf files (here: Loading *.adf files into QGIS?) was very good but I would like to know how to perform the same action using Python.


texts: EasyButton->how to change its location on map

texts:
I have a very simple question. I'm using EasyButton plugin of leaflet. But I can not change its location on map. For example how can I move the button to top right corner?

Here is my code:

   var toggle = L.easyButton ({
states: [{
stateName: 'remove-legend',
icon: 'masquer la légende',
title: 'masquer la légende',
onClick: function(control) {
    map.removeControl(legend0);
}
}, {
icon:'voir la légende',
stateName: 'add-legend',
  title: 'voir la légende' ,
  onClick: function(control) {
  if(map.getZoom()<12){
        map.addControl(legend12);
        }
}]
});
toggle.addTo(map);


texts: OpenGeo Suite create map application GXP or GeoEXT

texts:
I have seen a similar question but it never got answered. I am looking to develop a GIS system using the OpenGeo Suite, luckily I'm going to a conference on it in a few weeks, however that's weeks away.

I need to work on my front end, and I am really confused over GXP and GeoEXT, what should I be using to create my map application? the OpenGeo website has tutorials for both...

The system will enable panning/zoom, layer tree, maybe graphs, edit/features (click and change fields), grideditor, legend. Do I go down the OpenLayers route with GeoEXT or use SDK GXP. I am confused because some of the OpenGeo videos by themselves, seem to make out GXP is the way forward but then many examples and documentation they have, uses OpenLayers/GeoEXT.

This was the question, he basically wanted the same as me link


texts: Geoserver+postgis: how to complete creating database via REST API?

texts:
I'm using geoserver 2.13.1 with postgres 9.5 and postgis 2.4.2.

I am trying to use REST api to create workspace (works) and in it a datastore of postgis type. I have already prepared a template database in postgres and trying to use it here (this is the xml body of the post request to /geoserver/rest/workspaces/workspace_name/datastores):

<dataStore>
  <name>datastore_name</name>
  <description>dummy description</description>
  <type>PostGIS</type>
  <enabled>true</enabled>
  <connectionParameters>
    <entry key="Estimated extends">false</entry>
    <entry key="schema">venue</entry>
    <entry key="fetch size">1000</entry>
    <entry key="Max open prepared statements">50</entry>
    <entry key="Connection timeout">20</entry>
    <entry key="preparedStatements">false</entry>
    <entry key="database">mydbname</entry>
    <entry key="port">5432</entry>
    <entry key="min connections">1</entry>
    <entry key="dbtype">postgis</entry>
    <entry key="Loose bbox">true</entry>
    <entry key="host">localhost</entry>
    <entry key="max connections">10</entry>
    <entry key="user">foo</entry>
    <entry key="passwd">bar</entry>
    <entry key="create database">true</entry>
    <entry key="create database params">WITH TEMPLATE=mydbtemplate</entry>
  </connectionParameters>
  <__default>true</__default>
</dataStore>

The geoserver datastore is created as expected (and I'm getting 201 http code back), but the database itself is not created. No error messages.
Only if I open that store now via geoserver web interface and click on the save button (without changing anything), only then is the database really created using desired template. So, all those params are correct, but there is something missing on the server side, like a commit or something, that would trigger the create db job.

I have tried:

sending /reset and /reload requests afterwards - no db created
changing workspace and datastore properties with put request to trigger datastore saving somehow - changes are accepted, but no db created
I even tried sql injection by putting commit in that create database params entry after "with template" statement - also no error and no db created

I don't see any other REST API function looking even remotely related to "apply changes". So, how can I achieve that via REST API?


texts: QGIS Atlas web layer rendering - how can I slow it down?

texts:
I cannot get Atlas to export maps with web layers correctly.

When I look at the preview within composer everything works fine.
When I export an atlas the first page renders ok but for subsequent pages only a couple of tiles are usually generated.

My suspicion is that QGIS is rendering faster than the web layers can be downloaded.

I suggest that if it is possible to slow down the rendering process this will allow time for tiles to download.

I have no idea how to do this but maybe a developer can tell me how to slow down the Atlas rendering process to allow for tile downloads?

A similar question has been raised before on the GIS SE as Problems while rendering with Atlas.

My question is not the same because I am looking for a way to either test out my hypothesis or for somebody to tell me that my hypothesis completely wrong.


texts: In mapserver, how would I style polygons differently based on a separated points layer?

texts:
I have 2 tables in postgis - one of type polygons and one of type points.

I want the color of the polygons to be related to the number of points within each polygon, say white when no points are inside a polygon, pink when 'some' points are inside and red when 'a lot' of points are inside (a full range of colors based on percentage of points/polygon is what I'm really after).

I am using mapserver, and in the expression section of CLASS it doesn't mention any way of taking information from another table. Only attributes from within the same table used for the DATA tag.


texts: After having produced aspect of surface in GRASS GIS and being aware that GRASS GIS calculates aspect counterclockwise from east, how could we transform aspect from being a circular variable (0-360 degrees) into being a linear variable, so it can be used as explanatory variable in multiple regression models. I am aware of this answer (Transforming aspect values in order to use them in multiple regression models) which works perfectly in ArcGIS. How could the same procedure be carried out in GRASS GIS?


texts: Finding Lost Samsung S4 with Lat/Long Coordinates

texts: I lost my Verizon provider Samsung S4 phone while running on a narrow trail surrounded by brush.  The phone no longer has power.  Using Google Location History I obtained that Latitude and Longitude of the phone and have returned 2 times to the location looking for it using hard copies of assorted Google Map views and Bing aerial views centered on the coordinates.  I've yet to locate it in part due to the heavy cover.  Any suggestions?  My wife has an older Android phone powered by Virgin Mobile using Sprint's network.  My guess is if I brought that phone to the same coordinates it may give me a considerably different position (by up to 20 meters?) than the other phone.  Any thoughts on how I can find the phone?  Cell tower reception is also poor for both phones at the location where my phone was lost.


texts: Dynamic Segmentation Error for segments in Geomedia Pro 2015

texts:
I have a couple of segments that I am trying to do Dynamic Segmentation on with Geomedia Professional 2015.

I am using the EXACT same Route IDs for the segments and the EXACT same Postmiles for beginning and end segments. I am receiving the error message:

  Unable to locate event.  No matching key or wrong type fields are in
  the LRS dataset.

For one of the segments I created a new spreadsheet and even though the Route name is incorrect I got the segment to DynSeg.   I did that for the other two segments but it didn't work.

Will send segment database and jpg image of problems upon request.


texts: Drawing a Satellites Coverage FootPrint on an OpenLayers Vector Layer

texts:
We have been discussing offering at some point in the future satellite footprints on some of the maps we offer.  Initially the ISS but some Satellites may be appropriate as well. We have also been looking at this in the context of disaster management, and providing some tools for disaster planners to visualize where different useful satellites may be located to assist.

These two A and B excellent questions provide a great deal of information but are way too complex for my old brain to put into practice.

Would someone be able to show me the correct Openlayers JScript Code to draw the shape of the footprint on a vector layer, given the Satellites position and height.

Accuracy would not be critical as it is a visualization. Maps would be EPSG:4326 and/or EPSG:900913


texts: Alternative to Tippecanoe for generating Mapbox vector tiles on command-line

texts:
Tippecanoe is the tool I always use for converting Shapefiles/GeoJSONs to Mapbox vector tiles in PBF format, for rendering in Mapbox-GL-JS. But it doesn't seem ideal: its primary use case seems to be processing enormous point datasets such as tweets which means:

It focuses on points, not lines or polygons.
Many of its settings are focused on determining which features to drop.

There are other limitations in practice:

its processing of line and polygon features can be really ugly,
its "base zoom" model is not well explained, so that working through its dozens of configuration parameters is pretty hit and miss.
although it has some per-feature properties like maxzoom and minzoom it doesn't have a way to vary which properties are included per level
its functioning is pretty opaque, there's not much information generated about what choices were made at each level, how "heavy" the tiles are overall, etc.

So my question is: are there other ways to convert vector files to vector tiles that work better for polygons and lines?


texts:
I'm working on a project involving one database , but the customer's needs have evolved, so he needs to manage multiple projects. I tried to work with adding an id for each project, which gave me more complex queries all over the application, and I need to perform a replicate  per project , but the the replicates would be heavy for a one database containing multiple projects.

In this application I'm using multiple arcpy geoprocessing scripts, so if I will opt for multiple databases the SDE database would be a parameter.

But multiple databases would be easier to backup and to archive.

I tried to make a picture about the context of the question:

Question:

Is there any risk in adopting multiple geodatabases in one application?

Any advice or experiences are welcome.


texts: Extracting NDVI using QGIS from multiple images to Excel?

texts:
I created 25 images in QGIS.

Each image has different information about the NDVI of the area.

How can I extract mean, max, and min information about the NDVI from each image and put it in Excel?


texts: OGR and DXF files: Is it possible to handle layer names? (only "entities")

texts:
I'm using Python and OGR to extract the layers of a DXF file and convert them into SHP.

I began with:

import ogr, os, sys
driver = ogr.GetDriverByName('DXF')
datasource = driver.Open('test1.dxf', 0)
numberLayers = datasource.GetLayerCount()
for i in range(0, numberLayers):
    layer = datasource.GetLayerByIndex(i)
    layerName = layer.GetName()
    numberFeatures = layer.GetFeatureCount()
    print 'Layer=%s|Features=%s' % (layerName, numberFeatures)

but I get one only layer: entities (as expected: http://www.gdal.org/ogr/drv_dxf.html)

Obviously, it's the same information as if I execute: ogrinfo -so test1.dxf entities

My DXF file contains several CAD layers: "ROADS", "ELEVATION", ... Is it possible to handle these layer names with OGR? Thank you very much.


texts: Automatically write attribute when creating features in QGIS

texts:
Using QGIS 3.0 I have to draw a lot of features on the map and these features have attributes that change with each entry. However, there is a part of the attribute that is the same for all features. How do I automatically write this part of the attribute in the attribute table when creating a new point?

Example:

S1_00001
S1_00002
...
S1_0000N

The S1_ is always the same and can be in a different column as the serial numbering.


texts: Deleting file geodatabase with ModelBuilder when it has features in it gives ERROR 000601

texts:
I built a simple tool that creates 2 file GDBs for my work needs, and another to delete them after the job is done.

When the GDBs are empty, the delete tool works fine. However, when there are files in them, I can't use this tool to delete the GDBs, and this error pops:

If I manually delete all the files, then the tool works.

Suggestions?


texts: Extract raster values within shapefile with pygeoprocessing or gdal

texts:
I would like to know how to get all the raster values within a polygon using gdal or pygeoprocessing, without reading the entire grid as an array.

pygeoprocessing and gdal can do zonal statistics but only the min, max, mean, stdev or count are available from such a function. Since zonal statistics need to access the values, would it be easy to extract values the same way ?

I found a very similar question here :
(Getting pixel value of GDAL raster under OGR point without NumPy?)
but only for a particular "point".


texts: Extract the cells of a raster based on a logical query R

texts:
I have a raster Brick which represents the distribution models of 7 palm species, named currentStack_mask which looks like this.

As you can see I have 7 species and each of their rasters are represented with 0 and 1 values.

Basically what I want to do is (for each species) to extract all the cells that have a value of 1 and create another raster with those cells and of course to do it in R because I want to keep a track of what I am doing and also is because is faster and don't have to deal with all the intermediate files.

The equivalent function in Arcgis of what I want to do is Spatial Analyst Tools -> Extraction -> Extract by Attributes, which basically extract the cells of a raster based on a logical query, which in this case is that the cell value is 1.

I have tried with extract() function of the Raster package but this function extract the values not the cells.

Can anybody help me?... I am sure there is a short way to do this.


texts: Finding the shortest path with online services in QGIS

texts:
I am looking for a plugin that can detriment the shortest path between two points in QGIS.
It should not be a street line between the points, but the shortest way with consideration to the roads of the area in interest.


texts: Converting waveform LiDAR to discrete return format, in R?


texts: QGIS: How to exclude outliers from the legend classfication using an expression

texts:
I am trying to make an equal interval legend using a rate of burglaries per 1000 habitants using a combination of a column with the counts of burglaries and a column of total population using census districts. Due to the fact that there are outliers in the data (extremely high values) I would like to exclude them from the process of defining classes in the legend based on a threshold of the total population (less than 700).  I already build a expression in the expression function (which exclude values <700) window like this:

CASE WHEN  "burglaries_CS_POPN_2000|2000"  >= 700 THEN ( "burglaries_P500_ALL|2009"  /   "burglaries_CS_POPN_2000|2000" )  END

It works fine, but the polygons excluded are not showed anymore.  In ArcGis there is an option for "show symbol for excluded data" in the symbology-classification window.
Is there an option in QGIS that could perform the same task without adding a new column in the attribute table? Could it be done using an expression like I show above?


texts: How to get a Line string in OpenLayers 3?

texts:
I am using the code to get line string feature.when ever i am using coord2 array i am getting correct line string but whenever i am using points array i am not getting correct line string .I have observed that coordinate values in points array are stored in same format as in coord2 array.

This is my code:

   for(i=0;i<pathArr.length;i++)
      {
        var temp=(pathArr[i]+"").split(",");
      //console.log("lat lon "+temp[1]+"   " +temp[0]);
          points[i]= [temp[1],temp[0]];
     }
      //window.alert(points);

   var coord2=[[76,22],[77,24],[78,25],[66,12],[65,15],[76,22],[77,24],[78,25],[66,12],[65,15] ,[76,22],[77,24],[78,25],[66,12],[65,15],[78,25], [66,12],[65,15],[76,22]];
 window.alert(coord2);
  feature = new ol.Feature({
 geometry: new ol.geom.LineString(points)
     });


texts: Removing numeric characters from alphanumeric value in field calculation?

texts:
I have a field containing alphanumeric values in the following format:

A-AA00

Sometimes they are also

A-AAA0

I need a Python script I can put in the field calculator in ArcMap 10 that will at least remove the numeric characters, and retrain the alphabetical ones and the hyphen. So that, for example, this:

A-AA0

becomes this:

A-AA

I also want to replace the numeric chars with numeric chars from another field. This Stack Exchange entry addresses this issue, but the string has a space before the number and this can be used for splitting and indexing... my strings are not so handily constructed.


texts: CartoDB change public map in real-time with SQL query

texts:
I'm a data science student working on my capstone project. I'm new to CartoDB and have been playing around with their API a bit but can't seem to find what I'm looking for.

I have created a public map that I will display on my site, but I would like to have the map dynamically change based on a SQL query I send it.

For example, the dataset will contain the entire set of points, but I don't want all of them displayed at once, only I subset that I specify and change (through a select statement) as the user interacts with the program.


texts: Converting *.jpg scanned from line paper map to vector data in ArcGIS for Desktop?

texts:
I want to digitize automatically old paper map that scanned to jpg  .
I know that it's possible to digitize it with arc-map but I can not remember how to do it.


texts:
I am trying to run a buffer around a polyline (geography).  My first attempt below (based on ST_Buffer documentation:

SELECT ST_Buffer(ST_Geomfromtext(i30.geog),10,'endcap=round join=round');

gives an error:

missing FROM-clause entry for table "i30"
ST_Geomfromtext(**^here**i30.geog)


texts: How to edit a group of selected features in QGIS?

texts:
I am dealing with a large number of planning units across a very large study area.  I am selecting a number of features using the advanced search in the attribute table, and would like to edit them (i.e. change their colour).  The only way I can think of doing it is to save the selected features as a new layer, and then edit that.  Is there an easier way?
Any advice would be appreciated!


texts:
I'm using ArcGIS 10.2 and trying to obtain values in a field using SearchCursor.  Would someone help me with this?

import arcpy

oc = r"C:\Users\Scratch\oc.shp"

catchments = arcpy.da.SearchCursor(oc,"FEATUREID")

catchment = catchments.next()

for catchment in catchments:

   print catchment[0]

   print catchment[1]

catchments.reset()

The catchment[0] works, but if I try anything higher, I get tuple out of range error.
There are 11 rows in the shapefile, all with a unique FEATUREID.  I can only extract either the 1st row or all of them at once.  I need to be able to specify througout the rest of my code when I want the 1st, 2nd, 3rd, or 10th value in this shapefile under field = FEATUREID.


texts:
QGIS 2.18.2 from http://www.kyngchaos.com/software/qgis. Just reinstalled all packages, GDAL complete, NumPy, and matplotlib. Then reinstalled the app.

I created a polygon shapefile by drawing. But I don't see under Raster any options to clip. I want to clip a georeferenced scanned map.

http://www.qgistutorials.com/en/docs/raster_mosaicing_and_clipping.html mentions Raster ‣ Extraction ‣ Clipper. I don't have this option.

Other tutorials show similar options.

My Raster menu only have three options.

What am I missing?


texts: How to debug "Unable to load GDAL library" error in QGIS 2.0?

texts:
I have a problem with r.surf.contour. I create a DTM with v.to.rast.attribute and then use the r.surf.contour and it works fine in QGIS 1.9 but with QGIS 2.0 doesn't work.

I have a contour shapefile with elevation values. I use the v.to.rast.attr to rasterize the contour, but some pixels have null values, so I try to use the r.surf.contour to obtain a surface but in the end I've got an error. This is an excerpt: GRASS_INFO_ERROR(3016,1): Unable to load GDAL library GRASS_INFO_END(3016,1) C:\PROGRA~1\QGISDU~1\bin>g.region rast=output91bad3077faa4aa784bcedb42867ff89 And then GRASS_INFO_ERROR(1012,1): Raster map not found.

I think the problem is in raster from grass. Because I try to create a slope map with r.slope.aspect and it give me the same error. Please could you help me? My windows is Windows vista.


texts: Incorporating NDVI and texture metrics in OBIA-Supervised Image Classification

texts:
Anyone knows step-by-step tutorials for incorporating NDVI and texture metrics in land cover classification with ArcMap?

I've been doing supervised classification of Sentinel-2 with the assistance of segments resulted from Object Based Image Analysis to help me select the training sample.

However, many sources suggest to incorporate NDVI in the analysis as well as texture metrics since the spatial resolution of the image is quite high. For this matter, i would like to know how do you 'incorporate' NDVI and texture metrics practically in the software. Are there any calculations? or is it only a matter of visual assistance? Can i do it in ArcMap?


texts: Removing the existing context menus from legendlayeractions QGIS 2.6.1 using python

texts:
a = QAction( u"My Action", iface1.legendInterface() )
iface1.legendInterface().addLegendLayerAction( a, u"My Plugin Menu", u"id1", QgsMapLayer.VectorLayer, True )   #adds action
iface1.legendInterface().removeLegendLayerAction(a) #removes action

removes user made actions . I want to take out all the inbuilt context menu actions


texts: copy feature class and retain selected fields with arcpy

texts:
I've got a featureclass (fc) that I want to copy, but only retaining a selected number of fields, let's say field 9, 11, and 12 from the total of 15 fields. I want to use arcpy and suspect that fieldmappings are somehow required for this. I currently have:

fm= arcpy.FieldMappings()
fm.addtable(fc)
sel= range(0,9)+[10]+range(13,16)
for i in sel:
    fm.removeFieldMap(i)

However, this gives:

RuntimeError: FieldMappings: Error in removing a field map from field mappings for RemoveFieldMap

I can't quite see what I'm doing wrong? Perhaps it would be nicer to somehow select the fields I want and add these to a new fieldMapping object? All suggestions are welcome, thanks in advance!

Update. As suggested below, I now have the following:

inFields= ["fA", "fB", "fC"]
inTab= 'testFc'
fm= arcpy.fieldMappings()
fieldMap= arcpy.FieldMap()
for field in inFields:
    fieldMap.addInputField(inTab, field)
fm.addFieldMap(fieldMap)
arcpy.featureClassToFeatureClass_conversion(testFc, outDir, testOutFc, field_mapping= fm)

This works, a new featureclass is created. However, it contains only field fA. The python output reads:

Executing: FeatureClassToFeatureClass L:/GEOProjects/2Working.gdb/testFc :/GEOProjecten outDir testOutFc
# "fA "fA" true true false 8 Double 0 0 ,First,#,
L:/GEOProjects/2Working.gdb/testFc,fA,-1,-1,
L:/GEOProjects/2Working.gdb/testFc,fB,-1,-1,
L:/GEOProjects/2Working.gdb/testFc,fC,-1,-1" #
Start Time: Thu Oct 24 10:21:42 2013
Succeeded at Thu Oct 24 10:21:58 2013 (Elapsed Time: 16.00 seconds)

It seems that ony fA is being copied to the new featureclass. Any suggestion as to why this can be?


texts: Calculating volume with Cut Fill tool in ArcGIS Desktop?

texts:
This might be a very simple question - but it is quite vital and I want to make sure I am thinking this through right.

I have used the cut/fill tool to calculate the difference in volume of a volcano cone. I have an attribute table, is it just a case of adding up all the -ve values to find the overall volume? (negative because the volcano has inflated i.e. I should be gaining volume).

So, This is an example of my attribute table. However, I think that I actually need both the negative and the positive values. So, are you saying that I would have to do (2*-110676.30...)+(1*17935.13...)+(1*-55080.16...)+etc?


texts: How to compute weighted centroids in PostGIS?

texts: I have a group of points and each point has an associated count value. How can I  calculate the centroid of the group weighted by the values associated to each point?


texts:
I'm still sort of new to openlayers, but I've been able to set up a map which reads in a KML, loops through the different regions contained therein, and draws those polygons on a vector layer overtop of an OSM layer.  This all works fine until I try to add a control.  All I want to do is trigger an function when the polygon is clicked.  For testing purposes, right now I'm just trying to trigger an alert.
The applicable pieces of my code are below:

function onFeatureSelect(feature){
    alert("boom");
}

function onFeatureUnselect(feature){
    alert("unBoom");
}
var selectControl;
map = new OpenLayers.Map(divName, {
    controls: [new OpenLayers.Control.Navigation(), new OpenLayers.Control.LayerSwitcher(), new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.Attribution()],
    projection: googProj,
    units: 'm',
    maxResolution: 156543.0399,
    numZoomLevels: 10
});
var layerOSM = new OpenLayers.Layer.OSM();
map.addLayer(layerOSM);
var layerVector = new OpenLayers.Layer.Vector("KML Regions");

$.each(rPoly, function(){  //Looping through the data
    var feature_polygon = new OpenLayers.Feature.Vector( //We'll make a polygon from a linear ring object, which consists of points
new OpenLayers.Geometry.Polygon(new OpenLayers.Geometry.LinearRing(this)), {
        "rName": thisR
    }, {
        'fillColor': thisC,
        'strokeWidth': borderWidth,
        'strokeColor': borderColor,
        'fillOpacity': thisO
    });
    layerVector.addFeatures([feature_polygon]);
});

map.addLayer(layerVector);
selectControl = new OpenLayers.Control.SelectFeature(layerVector, {
    onSelect: onFeatureSelect,
    onUnselect: onFeatureUnselect
});
selectControl.activate();
map.addControl(selectControl);

Please remember that everything renders fine, I just can't get a function to fire when a polygon is clicked.
Any help would be GREATLY appreciated.  I'm getting pretty frustrated since this should be something fairly easy to do.
Thanks.

UPDATE

When I try using getLeyersByName to get the "KML Regions" layer, it comes back empty. How is this possible since all of the polygons are visible?  Could this be part of the issue?

UPDATE OF SHAME

What it actually turned out to be is that I had the loop creating the features on the layer running from a function and then I'd pass the fully rendered layer back. This was causing an error. So, in trying to simplify the code for the question, it seems like I fixed it if I would just paste it back in.
The code would have died anyway due to the calls being out of sorts, so I marked that answer as solving the problem.  Thanks everyone.


texts: Solve wrong topological error after union?

texts:
I have a shapefile which has accumulated quite a few topological errors (>400). Most of them are gaps and overlaps, which I know how to solve. However, when I use the "Topology Checker" Plugin checking my layer for gaps, it also identifies the entire extent of my polygons as being a gap (see image).

This is clearly a mistake and I am worried that these problems in the shapefile are the reason for me not being able to use other processing tools (e.g. the difference tool does not work properly on this layer).

How can I fix or understand in detail what the reason for this mistake is?

I am beginner at GIS and I am particularly interested to find out where I am using a wrong mindset or why am I creating these errors.

This is what I have done until now:
I have had an original shapefile without this error. I used union in ArcView 3.3. (because using Union in QGIS at that point led to many ruptures in the polygons and splinter polygons).

Is there a problem in compatibility in moving shapefiles between the two programs?

EDIT:
I have been trying to get around this error and started again from the original data with only one invalid geometry. After using union many errors appeared again (but not the one in the image). Then I used v.clean from Grass and this solved many other problems, but the error from the image showing a gap for the entire area reappeared. I am confused...


texts: Error in Mapproxy Tile generation for some map zoom level

texts:
I am using mapproxy to generate my map tile with wms configration as source, which is working fine for certain zoom level. But for zoom level like from 3 to 8 the map tile is not getting fully generated even waiting for long time.
Giving the tile output here

To get rid of this problem, I thought of implementing mapproxy-seeds configration, But getting some error like this:

Start seeding process (1 task)
  myseed1:
    Seeding cache 'latlong_cache' with grid 'webmercator' in EPSG:3857
    Limited to coverage in: -180.00000, -85.05113, 180.00000, 85.05113 (EPSG:4326)
    Levels: [3, 4, 5, 6, 7, 8]
    Overwriting: no tiles
[15:53:22]  3   0.00% -20037508.34279, -20037508.34279, 20037508.34279, 20037508.34279 (0 tiles)

An error occured. Retry in 2 seconds: IOError(13, 'Permission denied'). Retries left: 100

An error occured. Retry in 2 seconds: IOError(13, 'Permission denied'). Retries left: 100

An error occured. Retry in 4 seconds: IOError(13, 'Permission denied'). Retries left: 99

An error occured. Retry in 4 seconds: IOError(13, 'Permission denied'). Retries left: 99

An error occured. Retry in 8 seconds: IOError(13, 'Permission denied'). Retries left: 98

An error occured. Retry in 8 seconds: IOError(13, 'Permission denied'). Retries left: 98

An error occured. Retry in 16 seconds: IOError(13, 'Permission denied'). Retries left: 97

An error occured. Retry in 16 seconds: IOError(13, 'Permission denied'). Retries left: 97

An error occured. Retry in 32 seconds: IOError(13, 'Permission denied'). Retries left: 96

Here I am giving my seed.yaml and mapproxy.yaml configration.
seed.yaml

seeds:
  myseed1:
    caches: [latlong_cache]
    levels:
      from: 3
      to: 8

mapproxy.yaml

`services:
  demo:
  tms:
    use_grid_names: true
    origin: 'nw'
  kml:
      use_grid_names: true
  wmts:
  wms:
    md:
      title: MapProxy LatlongWMS Proxy
      abstract: This is a minimal MapProxy example.
    srs: ['EPSG:4326', 'EPSG:3857', 'EPSG:31467']
    bbox_srs: ['EPSG:4326', 'EPSG:3857', 'EPSG:31467']

layers:
  - name: latlong
    title: Latlong WMS - latlong.in
    sources: [latlong_cache]

caches:
  latlong_cache:
    grids: [webmercator]
    sources: [latlong_wms]

sources:
  latlong_wms:
    type: wms
    req:
      url: http://ip_address:81/cgi-bin/moon_map/qgis_mapserv.fcgi
      layers: world,landcovers,area,states,intersections,roads,way_points
      transparent: true
      transparent_color: '#d4d4c8'

grids:
    webmercator:
        base: GLOBAL_WEBMERCATOR`

Please anyone suggest me any solution for this.

One other isuue is that initialy at any zoom level the mapproxy is generating tile at very slow rate and taking too much time to load.So how can i fasten the initial map tile load.

How can I get rid of all these issues?


texts:
I have a table called Stores with 10k rows which has a column location geometry(Point,4326).

CREATE INDEX "Stores_spatial_gix"
ON "Stores"
USING gist
(location);

Upon KNN query

explain analyze select *
 from "Stores"
 order by ST_distance_sphere(location,st_point(-82.373978, 29.633657)) limit 3

I get about 800ms each time. What am I doing wrong that it takes so long time?


texts: Fire Danger Forecast index- classification - USGS

texts:
I recently downloaded a raster image called Fire Danger Forecast index from USGS - Fire Danger Forecast, and I wanted to get the probability of fire for a set of points I have in another layer, but the thing is that the symbology table does not comes with the data.

I find this very strange because without the table, it is very difficult to classify the whole range of color pixels (0-255) to get the fire probability (0%-100%) for each one.

There is a table available at the website but is a PNG (image) so I can't know where to break exactly.

Does anyone has any workaround about this issue?


texts: Adding active content to the Atlas in composer $scale

texts:
I have an Atlas of 178 pages that centers on different size villages, its all well organised through the composer and all works well. However, I can't use the scale bar as the differences between village area is far too great, I thought it would be simple just to add some rich text a link to the $scale of the sheet in the atlas..

Any ideas...?


texts:
I have a very simple setup where I have a feature vector layer on top of an OSM layer in OpenLayers 2.13.1. However, I notice that the icons of my points in the vector layer begin to 'become invisible' despite their SVG DOM bounds being positioned where the icon should be, making them sort of invisible as I zoom in.

The layer that holds my features is not overly exciting:

var ooiLayer = new OpenLayers.Layer.Vector('OOI Layer');

And here is how I create my feature vectors:

// Declared elsewhere but added for completeness sake
var ooiStyle = $.extend({}, OpenLayers.Feature.Vector.style['default'], {
    graphicWidth: 25,
    graphicHeight: 25,
    fillOpacity:.7,
    externalGraphic: 'img/ooi.png'
});

function add(ooi) {
    var wkt = new OpenLayers.Format.WKT(ooi);

    // The following line looks more complicated than it is. Assume the variable holds something like "POINT (31.123 34.321)" after evaluation.
    var wktData = ooi.entityInstancesGeometry[0].geometry.geometry.wellKnownText;
    var vector = wkt.read(wktData);
    if (ooi.entityTypeId == 14 && vector.geometry instanceof OpenLayers.Geometry.Point)
        this.createArea.call(this, vector.geometry.x, vector.geometry.y, ooi);
    else {
        var actualVector = new OpenLayers.Feature.Vector(
            latLon(vector.geometry.x, vector.geometry.y), // just a helper function that applies the correct projection if required; works as intended, so no worries here
            ooi,
            $.extend({}, ooiStyle, {
                title: ooi.entityName,
                externalGraphic: graphicFor(ooi.entityTypeId) // returns an URL to a 25x25 PNG
             })
        );

        ooiLayer.addFeatures(actualVector);
    }
}

From the looks of it, it almost appears as if the first icon on the map sets the visual bounding box for all the other icons. I've created a small example to illustrate this behavior: http://youtu.be/_axFRI6Pw1U

What is wrong here?


texts: Calculating standard metrics in multiple LiDAR files with lidR 2.0

texts:
I am porting my existing code to lidR version 2.0 and would like someone to let me know if I am coding efficiently.  My objective is:

to process a directory of LAS files
normalize las files
calculate stdmetrics using grid_metrics
write out results in a data frame that I can use for subsequent analysis

To perform all these step at once I'm trying to use the function catalog_apply to apply user-defined function. What I did so far is:

Load a LAScatalog and set some processing options

las.catalogue <- catalog("F://Work//LASTest")
opt_cores(las.catalogue) <- 4L
opt_chunk_size(las.catalogue) <- 260
opt_chunk_buffer(las.catalogue) <- 10

Define my function that is expected to read sequentially some chunks of the point cloud, normalize this point cloud chunk and computed metrics in a area based approach.

stat.function <- function(las)
{
  las = lasnormalize(las, tin())
  metrics <- grid_metrics(las, .stdmetrics_z, res = 20)
  # I think I am supposed to remove buffers but need advice on how to do
  return(metrics)
}

Then with catalog_apply I sequentially run my function on my catalog

metrics <- catalog_apply(las.catalogue, stat.function, .options = opt)

To finish I merge the list of outputs for each chunk into a single RasterBrick. Because I want a data.frame I convert the RasterBrickinto a data.frame

metric.merge  <- do.call(raster::merge, metrics)
metrics <- raster::as.data.frame(metric.merge,xy=TRUE, na.rm = TRUE)

I'm not yet comfortable with the new design of the lidRpackage. Am I on the right track to perform my computation?


texts: Updating apps on OpenGeo suite Community Edition

texts: I would like to update PostgreSQL and PostGIS on OpenGeo suite Community Edition to a newer version, maybe the other applications too (GeoServer, OpenLayers, etc.). Can it be done? Can someone point me the right direction?


texts:
I am having trouble collecting GPS coordinates. In my case, I need to walk from point A to point B and collect the coordinates avoiding any obstacles. I used an android phone and mapmyrun app. When I start from point A, I can see it does not start from where I am standing. It shows I am at the other side of the road. While walking along the footpath sometimes the GPS takes me other side of the road and shows my location in a restaurant or in a building whereas I am walking in the footpath. Now I have the coordinates of these walks. How can I remove these errors from the coordinates?
Note I also used other apps like adidas Runtastic. It gives me same problem. The GPS location changes a lot.


texts: Convert Blob file to an ArcObjects object

texts:
I have a blob file which I take from a database. This blob is actually an MXD file. I want to load this MXD to an axLayoutControl (ArcGIS Engine). But for security reasons I do not want to write this on disk and then load it with LoadMXD(string path) method. I want to pass it to my layout control directly. So I need to convert it to some interface like IMXDContents. Below is the code I use. But the exception "Class not registered (Exception from HRESULT: 0x80040154 (REGDB_E_CLASSNOTREG))" is thrown on ObjectStream.LoadObject(ref guid, k) . I googled this and the only thing I found was that I should change my platform target to x86 instead of Any CPU. Believe me I tried this in every posible way. And still not working. Can anyone please help?

IXMLStream pXMLStream = new XMLStreamClass();
pXMLStream.LoadFromBytes(ref bytes); //bytes is the cast of blob file to bytes array

IObjectStream ObjectStream = new ObjectStreamClass();
ObjectStream.Stream = pXMLStream as IStream;

IMxdContents k = new PageLayoutClass();
Guid guid = typeof(IMxdContents).GUID;

ObjectStream.LoadObject(ref guid, k);


texts: Search and zoom json features in Leaflet

texts:
How can I add a search box to a map that allows users to search>zoom to a polygon name?

So, if there are 100 polygons, the user would be able to type the name of the polygon (i.e., Symantec), then the map would zoom to the matching county.


texts: How to create a KML that supports animated ground overlays?

texts:
I have 12 images that represent weather conditions at day 1 through day 12. I am using gdal2tiles to produce map tiles from each of these images. gdal2tiles also produces a KML file that can be used in Google Earth to display the tiles.

However, I want to create an animated layer out of all 12 images.
Google Earth seems to support animation in KML files (see https://developers.google.com/kml/documentation/time#example2). Are there any tools to create an animated ground overlay KML file such as this? Is gdal2tiles capable of it?

Thanks a lot.
- Mike


texts: How to handle multiple rasters (precipitation data) in QGIS?

texts:
Currently I am working at a hydrological model for Rwanda which needs monthly precipitation data averaged over the area of subcatchments.

The precipitation data are 144 TIF's from CHIRPS for timespan 2004-2015.
The subcatchments are shapefiles

The output has to be a .txt file with the monthly precipitation data per subcatchment.

Is there a way in QGIS to group the 144 raster files and clip and calculate them together?

I use QGIS 2.12.1 with GRASS 6.4.3, Python 2.7 and PostGreSQL 9.6. Most of my experience is in QGIS and GRASS, and I am learning PostGIS since two weeks.


texts: Cookie cut a polygon out of another polygon using JTS only

texts:
I need to do some surgery on polygons in mongodb
I have a polygon I want to cut from
            //this is a test case, in reality i am using mongodb to store geoJsonPolygons..so converting to JTS Geometry here.
                Geometry parent = GeoJsonPolygonUtil.toGeometry(geoJsonPolygonParent);

                //Cookie cut this polygon from parent.
                Geometry child = GeoJsonPolygonUtil.toGeometry(geoJsonPolygon);

                //create line segments from the polygon we want to cut out..
                Geometry lineSegements = PolygonUtil.polygonize(child);

                //split?(cut) out the polygon..
                Geometry parentAfterChildRemoved = PolygonUtil.splitPolygon(parent, lineSegements);

                //create a new polygon minus the child to cut out...
                geoJsonPolygonParent = GeoJsonConversionUtil.extract(parentAfterChildRemoved.getCoordinates());

I am using the PolygonUtil polygonize, splitPolygon from another gis stack overflow answer thinking it would work...but it does not, the parent polygon has no changes...
JTS: split arbitrary polygon by a line
Parent Polygon - zipcode 92683
https://gist.github.com/boundaries-io/4769bc85338e7398dfdf27307b3d381d

Child Polygon - zipcode 92655
https://gist.github.com/boundaries-io/139e7645a59ef7a36aecfe7d42f9beac

my PolygonUtil code:
public class PolygonUtil {
    public static Geometry polygonize(Geometry geometry) {
        List lines = LineStringExtracter.getLines(geometry);
        Polygonizer polygonizer = new Polygonizer();
        polygonizer.add(lines);
        Collection polys = polygonizer.getPolygons();
        Polygon[] polyArray = GeometryFactory.toPolygonArray(polys);
        return geometry.getFactory().createGeometryCollection(polyArray);
    }

    public static Geometry splitPolygon(Geometry poly, Geometry line) {
        Geometry nodedLinework = poly.getBoundary().union(line);
        Geometry polys = polygonize(nodedLinework);

        // Only keep polygons which are inside the input
        List output = new ArrayList();
        for (int i = 0; i < polys.getNumGeometries(); i++) {
            Polygon candpoly = (Polygon) polys.getGeometryN(i);
            if (poly.contains(candpoly.getInteriorPoint())) {
                output.add(candpoly);
            }
        }
        return poly.getFactory().createGeometryCollection(GeometryFactory.toGeometryArray(output));
    }
}


texts: Error while using raster calculator in model builder

texts:
I am trying to use the model builder for the RUSLE Model. Well this is how the model looks like:

So in the raster calculator I have a formula which also involves the rasters generated in the model. This is error which I get .Could anyone tell me what could be going wrong.

The  formula which I am using in the first raster calculator block is :

Power( ["%FlowAcc_Flow%"] * 10 / 22.1, 0.4) * Power(Sin [" %Slope_Fill%" ]  / 0.09, 1.4) * 1.4


texts: ArcObjects 10.1 .NET Convert Point from WGS1984 to OSGB1936

texts:
Please can someone help me convert a point in WGS1984 to OSGB1936, I have followed all the documentation (or so i think...) and yet have not made any headway into this problem. This is my code :

Type factoryType = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
System.Object obj = Activator.CreateInstance(factoryType);
ISpatialReferenceFactory3 spacialReferenceFactory3 = obj as ISpatialReferenceFactory3;

// Create Transformation from WGS84 to OSGB86
IGeoTransformation geoTrans = spacialReferenceFactory3.CreateGeoTransformation((int)esriSRGeoTransformationType.esriSRGeoTransformation_OSGB1936_To_WGS1984Petrol) as IGeoTransformation;

ISpatialReference fromSpatialReference;
ISpatialReference toSpatialReference;
geoTrans.GetSpatialReferences(out fromSpatialReference, out toSpatialReference);

IGeometry5 geometry;
IPoint point = new PointClass();
point.PutCoords(51.465615, -3.159875);
geometry = point as IGeometry5;
geometry.SpatialReference = toSpatialReference;
geometry.ProjectEx(fromSpatialReference, esriTransformDirection.esriTransformForward, geoTrans, false, 0.0, 0.0);
point = new PointClass();
point = geometry as IPoint;

My point is being returned from the Geometry object as 51.46200954514034, -3.1549783406850165. I have tried using the reverse transformation as specified in geometry.ProjectEx(), and all i get is the same X and Y as the original point. I have also tried different combinations of reverse / forward and also swapping the fromSpatialReference and toSpatialReference.


texts: QGIS geoprocessing scripts, default values for field names

texts:
In QGIS geoprocessing scripts, default values must be provided for input numbers and strings, according to http://docs.qgis.org/2.8/en/docs/user_manual/processing/console.html#creating-scripts-and-running-them-from-the-toolbox like:

##mynumer=number 3.1415329

Since I often have to process layers with all the same field names (unfortunately not assurable in any case), I would like to pass a default value for a field name to the geoprocessing script, following a logic like 'take the default field name if it exists in the layer referred to, otherwise let the user select one'.

Is this possible somehow, perhaps with a bit of python?


texts: How to speed up 'intersection' tool for big amount of data?

texts:
I can't figure out how to work with a large dataset. The problem is: I want to know how many points (over 700 000 adresses) are within a polygon layer (it's made out of 5 layers, train noise)(they are both shapefiles). Using the 'Join by location' tool or 'Add polygon attributes to point' (SAGA) would take days to calculate it. I tried using a Spatial lite database and a Postgis database (SQL) but that also took way too long (10 attributes took over a minute to calculate). This is the script I used:

SELECT
bag_panden."id", trein_geluid_nacht."fromto"
FROM
bag_panden
inner JOIN
trein_geluid_nacht
ON ST_Within(bag_panden.geom, trein_geluid_nacht.geom)
where
bag_panden."id" <= 100

If anyone knows a way to calculate this faster I would like to know how. Maybe a different tool? Some Settings? A different script?

I am using QGIS 2.18.14


texts: Making QGIS layer update from changed .geojson data source

texts:
I have the same problem as Making QGIS layer update from changed data source.
However, I cannot get the solution proposed by Matthias Kuhn in a comment on an answer to that question to work.
I am using QGIS 3.4, I have loaded a .geojson file as a vector layer and am trying to refresh it without removing and re-adding. The map refreshes when I move it but the attribute table is not re-loaded.
I may not have understood how to use the reload function so here are my steps, I am using the Python Console.
       mc = iface.mapCanvas()
       layer = mc.currentLayer()
       layer.name()
       'radarOutGJTemp'  #--- Correct

       # Tried the following with errors attached
       layer.dataProvider().reload()
AttributeError: 'QgsVectorDataProvider' object has no attribute 'reload'

       layer.dataProvider.reload('radarOutGJTemp')
AttributeError: 'builtin_function_or_method' object has no attribute 'reload'

       layer.dataProvider.reload(radarOutGJTemp)
AttributeError: 'builtin_function_or_method' object has no attribute 'reload'

I can see the function reload appear in a window as I type. I did try to import reload just in case, that failed but again I could see it appear in the window as I typed.
import reload
ModuleNotFoundError: No module named 'reload'


texts: Problem when running the Clip Raster by Mask Layer - Layers were not correctly generated

texts:
I need help, this is for my Final Year Project.
Yesterday before I did something on my QGIS, I can run the mask layer but suddenly after I followed a tutorial on youtube on how to do 3D Map, I cannot run the mask layer again. I've search answers but none of them can help me. What should I do? Below is from the log messages.
QGIS version: 3.16.2-Hannover
QGIS code revision: f1660f9da5
Qt version: 5.11.2
GDAL version: 3.1.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'Clip raster by mask layer' starting…
Input parameters:
{ 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : 'C:/Users/User/Desktop/New folder/Peta/DEM/Malaysia.tif', 'KEEP_RESOLUTION' : False, 'MASK' : 'C:/Users/User/Desktop/New folder/Peta/Shapefile/Batu Gajah/Shale.shp|layername=Shale', 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : None, 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None }

GDAL command:
gdalwarp -of GTiff -cutline "C:/Users/User/Desktop/New folder/Peta/Shapefile/Batu Gajah/Shale.shp" -cl Shale -crop_to_cutline "C:/Users/User/Desktop/New folder/Peta/DEM/Malaysia.tif" C:/Users/User/AppData/Local/Temp/processing_FQruTL/38b9353213c7499cb36cab6f30c2a177/OUTPUT.tif
GDAL command output:
Warning 1: Self-intersection at or near point 100.9829299492848 4.3754446642990406

ERROR 1: Cutline polygon is invalid.

Execution completed in 0.16 seconds
Results:
{'OUTPUT': 'C:/Users/User/AppData/Local/Temp/processing_FQruTL/38b9353213c7499cb36cab6f30c2a177/OUTPUT.tif'}

Loading resulting layers
The following layers were not correctly generated.
• C:/Users/User/AppData/Local/Temp/processing_FQruTL/38b9353213c7499cb36cab6f30c2a177/OUTPUT.tif
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: How to make arcpy to iterate through folders & sub-folders to access gdb's?

texts:
With reference to Check whether any feature classes in multiple gdb's has feature class representation, the code I've mentioned in the other thread is working great for single folder alone. We have 100's of gdb's in different folders & sub-folders.

Since I am very much new to arcpy, I don't know how to do this. I searched & found arcpy.da.walk has to be used to iterate through folders. I've tried that (code below). The code is getting processed but no results or not even throwing any errors. But the process gets completed every time.

Edit:

Based on the comments and code given in this link, i've changed the code

import os
import arcpy
def FindField(fc,myField):
    fieldList = arcpy.ListFields(fc,myField)
    for field in fieldList:
        if str.lower(str(field.name)) == str.lower(myField):
            print gdb, fc + " contains fieldname: " + myField
myField = "RuleID"
top_folder = r"W:\RT_QAQC\Received\20160411_GDM_QA\VectorDelivery\RegionalReference"

for path, dirs, files in os.walk(top_folder):
    for d in dirs:
        if not d.endswith(".gdb"):
            continue
        gdb_path = os.path.join(path, d)
        arcpy.env.workspace = gdb_path

        for fc in arcpy.ListFeatureClasses():
            FindField(fc,myField)

        for fds in arcpy.ListDatasets('','feature'):
            for fc in arcpy.ListFeatureClasses('','',fds):
                FindField(fc,myField)


texts:
I am trying to support a WMS layer with a North Polar Stereographic projection (EPSG3995) on Leaflet but cannot make it to work.

I used Proj4Leaflet and have created the projection with the resolutions given by GeoServer (the WMS layer is served by GeoServer, see joined image).

On the following jsfiddle: https://jsfiddle.net/gaubert/Lyojcrzo/
I have tried to implement it and if you zooming on the leaflet map, you will see that the tiles are out place?

What is missing to make it work? Do I need to implement a transform or project method?


texts: Make user-defined feature set stops first and last in sequence using ArcGIS ModelBuilder and Network Analyst?

texts:
I have about 800 stops, in 26 category groups (A to Z). Each category group represents a 'run' or 'route'. What I want my model to do is allow the user to define (by clicking on the map) a start and end location and then solve one of the routes. So when the model is run...

The dialogue box should ask the user for the Run Name and then to click two locations on the map. Then, it should somehow make the two user-selected locations the 1st and last stops in the sequence and solve the route.

Any ideas?


texts:
Trying to use QGIS 2.0 on Mint 15 (based on Ubuntu 13)

I installed, uninstalled and re-installed QGIS from various repositories (ubuntugis-unstable, then qgis.org/debian raring - not clear for me which one I should use). I don't seem to have any python support. Any help appreciated!

If I try

apt-get install python-qgis

I get the message that

the packages libgdal1 libqgis2.0.1 qgis qgis-providers
are going to be uninstalled

Bump:
Still no solution. I uninstalled QGIS and apt-get autoremove, kept only ppa:ubuntugis/ubuntugis-unstable. Now when I try

sudo apt-get install qgis python-qgis qgis-plugin-grass

the error message is

Les paquets suivants contiennent des dépendances non satisfaites :
 python-qgis : Dépend: python-qgis-common (= 2.0.1-2+raring1) mais ne sera pas installé
 qgis-plugin-grass : Dépend: grass642
E: Impossible de corriger les problèmes, des paquets défectueux sont en mode « garder en l'état ».

Can someone help?


texts: Is there easy way to get layers names from WMS server?

texts:
I've been looking for information on how to get layers names from certain WMS servers, but I didn't find any.

E.g. I have this server: http://sdi.geoportal.gov.pl/wms_prg/wmservice.aspx and I want to know what layers can I get from it. Is there any chance to know their names?

Thanks in advance!


texts:
I've noticed there is not an option to make the cell size for both X and Y the same when interpolating point data in QGIS 3.0.

When adjusting the cell size in QGIS 2.18, the columns would automatically adjust. Someone once told me the more ocean that is in the frame, the more askew the resultant interpolated raster will be, unless the cell sizes are set the same.

Below is a screen shot from 2.18

And from 3.0

This is an important feature for me because if the resultant raster pixels end up rectangular instead of square, any viewshed analysis I perform end up some 15-50km to the north of my area of interest.


texts: What is the Scan Direction of a satelite image?

texts:
There are two possible values that I am aware of for describing the scan Direction of a satellite capture:

Forward
Reverse

What is the meaning of  Forward and Reverse scan direction?
What are the pros and cons of each direction?
Provided I have a have choice should I favor one from the other? Which one?


texts: Removing cloud cover on LandSat 8 imagery using ArcGIS for Desktop?

texts:
How can I remove cloud cover on landsat8 imagery?
using ArcGis10.2


texts: Centering Feature before exporting map to PDF?

texts:
I am trying to export single point file maps to PDF. My maps are exporting fine, but I'm having trouble centering the point file in the middle of the map, it's always outside of the data frame. Any suggestions?

import arcpy

mxd = arcpy.mapping.MapDocument(r'P:/Python/test.mxd')
substation = 'Mcarthur.shp'

add_substation = arcpy.mapping.Layer('mcarthur.lyr')

df = arcpy.mapping.ListDataFrames(mxd, 'Data_Frame')[0]

for df in arcpy.mapping.ListDataFrames(mxd):
  df.rotation = 0
  df.scale = 250000

sub = arcpy.mapping.AddLayer(df, add_substation, 'AUTO_ARRANGE')

north_arrow = arcpy.mapping.ListLayoutElements(mxd, 'MAPSURROUND_ELEMENT',
'North*')
narrow = north_arrow[0]
narrow.elementPositionX = 4.25
narrow.elementPositionY = 1.25

map_legend = arcpy.mapping.ListLayoutElements(mxd, 'LEGEND_ELEMENT',
'Map_L*')
legend = map_legend[0]
legend.elementPositionX = 0.5
legend.elementPositionY = 1.33

map_scale_bar = arcpy.mapping.ListLayoutElements(mxd, 'MAPSURROUND_ELEMENT',
'Scale*')
scalebar = map_scale_bar[0]
scalebar.elementPositionX = 4.25
scalebar.elementPositionY = 0.33

map_title = arcpy.mapping.ListLayoutElements(mxd, 'TEXT_ELEMENT', 'Map_T*')
title = map_title[0]
title.elementPositionX = 4.25
title.elementPositionY = 10.75
title.text = substation

map_data_frame = arcpy.mapping.ListLayoutElements(mxd, 'DATAFRAME_ELEMENT',
'Data_F*')
dataframe = map_data_frame[0]
dataframe.elementPositionX = 4.25
dataframe.elementPositionY = 2.75
dataframe.elementHeight = 7.5
dataframe.elementWidth  = 7.5

mpc_image = arcpy.mapping.ListLayoutElements(mxd, 'PICTURE_ELEMENT', 'MPC*')
mpcimage = mpc_image[0]
dataframe.elementpositionX = 8.0
dataframe.elementpositionY = 1.55

copy = mxd.saveACopy(r'P:/Python/mapoutput.mxd')
arcpy.mapping.ExportToPDF(mxd, r'P:/Python/mapoutput.pdf')

del mxd


texts: Turn off reading zip archives in browser

texts:
I've recently learned that QGIS Browser Panel is able to read zip archives.

Unfortunately when there's a large archive (let's say 10GB) in a folder the browser is forced to scan it's contents. This prevents display of the folder's content until the scan is finished.

Is there a way to turn off zip archive read by the Browser Panel?


texts: Create Several Feature Layers and add them to a MXD

texts:
Every time I run this script now, I get an error on the line with the MakeFeatureLayer line. Sorry for all the questions! I'm really new to Python.

import arcpy, os

arcpy.env.workspace = r'C:\pathway to gdb'
fcs = arcpy.ListFeatureClasses()

mxd = arcpy.mapping.MapDocument(r'C:\pathway to mxd')
DF = arcpy.mapping.ListDataFrames(mxd,'*')[0]

for feature in fcs:
    inputlayer = 'C:\pathway to gdb'
    outlayer = os.path.join(fcs, 'feature'+'.lyr')

    arcpy.MakeFeatureLayer_management(inputlayer,outlayer)

    layer = arcpy.mapping.Layer(mxd)
    arcpy.mapping.AddLayer(DF,layer,"AUTO_ARRANGE")


texts: Determine density or cover from a raster layer that has been clipped by a polygon

texts:
I am a very new GIS user using ArcMap 10. I have a raster layer of tree canopy that I have clipped so that the layer is in a polygon that is a circle with a 2km radius.

I have done this for four circles and I am trying to compare the tree canopy across the four circles. Something like a density or % area covered would be perfect.

The problem is... the tree canopy raster layer attribute table only has counts for where the canopy is... and no counts for where it is not. I could get a % cover (for comparisons across the four circles) if I knew what the count of the no tree canopy cells was, but this information is absent. It is almost as if the raster layer is like a point layer, but I cannot count each cell as a point. I hope this makes sense. I am completely lost.. any help would be appreciated!

Thanks in advance!


texts: Snapping two line shapefiles according to attribute field value using ArcGIS Desktop?

texts: I have a two line shapefiles (Trace.shp and Snap.shp). There is a GIS_ID field in both files.  Now I want to snap "Snap.shp" to "Trace.shp" according to GIS_ID value. I have tried to use Snap tool for ArcGIS Desktop but I am not sure about the results because using that solution I am not sure that line is snapped to the right line from another layer with the same GIS_ID. I want to do that using ArcGIS 10.4


texts: Algorithm fails (not found) : gdalogr:warpreproject

texts:
I am using the code in "QGIS Python Programming Cookbook", testing the code to learn about programming in QGIS. But when I try to reproject a raster, I get the error: Error: Wrong number of parameters. When I type in processing.alghelp("gdalogr:warpreproject") it returns "Algorithm not found".  Here is my code:

import processing
from processing.core.Processing import Processing
from qgis.core import *
from qgis.gui import *

Processing.initialize()
Processing.updateAlgsList()
#load a raster layer and set a COORDINATE SYSTEM
rasterLyr = QgsRasterLayer("G:\...\INPUT.ecw","MTN50")
print rasterLyr.isValid()
crs = QgsCoordinateReferenceSystem(25830, QgsCoordinateReferenceSystem.EpsgCrsId)
rasterLyr.setCrs(crs)

QgsMapLayerRegistry.instance().addMapLayers([rasterLyr])
#now, REPROJECT
processing.runalg("gdalogr:warpreproject", rasterLyr, "EPSG:25830", "EPSG:3722", None, 0, 0, 0, 0, 2, "\0, C:\..\output.ecw")

I am using QGIS Essen, 2.14. Is it a problem of using a code that isn't updated?


texts: Using arcgis.rand() with variables from code block of Calculate Value tool?

texts:
Within a larger ArcGIS ModelBuilder model I am trying to populate a field in an attribute table with a random number, using Calculate Value. The overall distribution of random numbers needs to follow a Normal distribution. The parameters of that normal distribution (mean and standard deviation) need to be calculated from model variables specified by the user. Is this possible?

The arcgis.rand() function can do what I want and works perfectly if I specify the mean and stdev directly in the Calculate Value expression:

arcgis.rand("Normal 6 2.5")

If I use model variables directly, it also works:

arcgis.rand("Normal %mean% %stdev%")

However, in some instances I need to do some simple calculations to obtain the parameters from variables the user has entered. For example, they have already entered variables A, B and C, and I need the mean to be A*B and the stdev to be B/C.

I have attempted to do this via the code block:

import random
MN = %A%*%B%
SD = %B%/%C%

and then the expression

arcgis.rand("Normal MN SD")

but this always results in values following the default distribution (mean 0, stdev 1). Troubleshooting shows this happens even if I define the code block variables as simple numbers, suggesting the problem seems to be using code block variables in the expression, not in the definition or calculation of the variables themselves.

My question is: how can I do what I need using arcgis.rand()?

I am open to suggestions of other methods/workflows for achieving what I need, but I am a novice at Python, in case that is not already obvious! I am using ArcGIS Desktop Advanced 10.2.2.


texts: How to batch-process attribute tables in QGIS (or GRASS)?

texts:
I'm using GRASS 6.4.3 in Linuxmint 16.
  I have a bash script for GRASS which I use to process almost 300 raster files.
The script first imports some raster files (r.in.gdal), then does some calculations and classification (using r.mapcalc and r.reclass). And finally transforms the resultant raster to vector (r.to.vect) and saves it (v.out.ogr, using type=area).

Let's say each vector files has a field labeled category and that it takes only two values 1 and 2.
As I'm interested only in the total area of polygons with category = 2 and with an area higher than 10,000 m2, I open each vector file in QGIS, edit its attribute table, calculate the area of each polygon and eliminate all of them with category=1 and $area<10,000.

But I have to do it manually,

Is there a way to do it in a batch process within QGIS?

Or even better, is it possible to do it in GRASS?

Is it possible to edit the attribute table within the command console in GRASS?

Could I do it in the same bash script which I use to do all the other processes in GRASS?

I've been reading this but I'm not sure how to apply it to my problem. I couldn't find how to calculate the polygons area and add it to a new column using GRASS.
I think I could use db.select to search for the polygons I want to eliminate, but I don't know how to do it.

Thanks for your help!


texts: Swift library to convert latitudes and longitudes from EPSG:4326 into EPSG:3857

texts:
I would like to know a Swift library or method that I can use to convert from 4326 to 3857. I have referred some libraries such as Proj4Swift and Proj.4, but to no avail. The first link does not has a workable download and the second does not has a documentation. Would anyone suggest a library or method that I can use to achieve this in Swift?

I have also referred other similar questions but could not find a proper solution.


texts:
I have a shapefile with many polygons in it. I want to dissolve by two fields.

Is there any possible way to dissolve this by NAMA_KPH and PL15_ID?


texts:
did anyone try rerendering OpenLayers in JSF2? I am using a4j:ajax events,
I mark the render atribute with render="map" where map is naturally the id
of the openlayers div tag, and I use oncompleate="init()" method to reload
the background data. And it works apart from the fact that the ''old' map
stays in place and the 'new' one is displayed belove. If I use the rerender
event again (i.e. changing some JSF components like selectOneMenu) than a
3rd map is added to the bottom, and so on and so on...

does any one know a
ways to fix this?


texts:
I need to georeference a list of places. Unfortunately, my original file doesn’t contain addresses, it is made up of a list of famous locations, such as the St Mark Basilica, Miramare Castle, Scrovegni Chapel.
Google my maps recognizes these famous places, and so, even if my file doesn’t contain addresses, Google My maps can localize these places correctly on the map (in Trieste, Venice, Padua) and creates markers.

When I export a KML or KMZ file from My Maps and open it on QGIS, I only get a table containing the information.

I was wondering if you knew any tools in QGIS that do that or any other tool in the net? I saw that there are several people asking this question, but my case is rather different as I don’t have any street address field.


texts: Performing Union_analysis on outputs created with arcpy.GetParameterAsText()?

texts:
I am trying to write a script tool for my geospatial programming class and I am stuck (this is an online course and we get NO support from our instructor).

What I need to do is create buffers (different values) around three different features, these are set up as parameters that the user will define. I got this to work just fine.

Next, I need to perform a Union_analysis on those three buffer outputs, this is where I am getting stuck.  I did this in a previous script, where I was using known feature classes, but I cannot figure out how to get this to recognize the outputs of the three "user selected" features.

Here is my script so far:

import arcpy

arcpy.env.overwriteOutput = True
buffered_area = "buffered_area"

# Get the input parameters for the Buffer tool
inPath = arcpy.GetParameterAsText(0)
outPath = arcpy.GetParameterAsText(1)
bufferDistance = arcpy.GetParameterAsText(2)
inPath1 = arcpy.GetParameterAsText(3)
outPath1 = arcpy.GetParameterAsText(4)
bufferDistance1 = arcpy.GetParameterAsText(5)
inPath2 = arcpy.GetParameterAsText(6)
outPath2 = arcpy.GetParameterAsText(7)
bufferDistance2 = arcpy.GetParameterAsText(8)

# Run the Buffer tool
buffer1 = arcpy.Buffer_analysis(inPath, outPath, bufferDistance)
buffer2 = arcpy.Buffer_analysis(inPath1, outPath1, bufferDistance1)
buffer3 = arcpy.Buffer_analysis(inPath2, outPath2, bufferDistance2)

arcpy.Union_analysis("buffer1 #; buffer2 #; buffer3 #", buffered_area, "ALL", "", "GAPS")

And here is the error message I get:

Failed script Script1...

Traceback (most recent call last):
  File "\\SNHU-LABS-UTIL2.snhu-labs.com\userprofiles$\william.herzig\Documents\programming\try.py", line 29, in <module>
    arcpy.Union_analysis("buffer1 #; buffer2 #; buffer3 #", buffered_area, "ALL", "", "GAPS")
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\analysis.py", line 629, in Union
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset buffer1 #;buffer2 #;buffer3 # does not exist or is not supported
WARNING 000725: Output Feature Class: Dataset \\SNHU-LABS-UTIL2.snhu-labs.com\userprofiles$\william.herzig\Documents\ArcGIS\Default.gdb\buffered_area already exists.
Failed to execute (Union).

Failed to execute (Script1).


texts:
I have a dataframe with x,y coordinates. I've sampled 10 x,y coordinates from my dataframe:

y <-c(24.637032,21.569813,21.687,19.677183,24.72438,21.433136,24.637032,22.21065,21.443762,25.273882)

x <-c(-156.02597,-112.00287,-112.00287,-105.27803,-159.47324,-157.858,-156.026,-109.68626,-112.12681,-112.15393)

I've converted these points to a SpatialPoints objects:

library(sp)
myPoints <- data.frame(x,y)
coordinates(myPoints) <- c('x', 'y')
proj4string(myPoints) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

I've used these points to make a SpatialGridobject:

bb <- bbox(myPoints)
myPointsGrid <- GridTopology(cellcentre.offset =  bb[, 1] + (cs/2), cellsize = c(1, 1), cells.dim = ceiling(diff(t(bb))/cs))
myPointsGrid <- SpatialGrid(myPointsGrid, proj4string = "+proj=longlat +ellps=WGS84 +datum=WGS84")

I need two things 1. find which grid cells contain a point 2. Find the number of points in each grid cell. I tried...

over(myPointsGrid, myPoints)

...but this has thrown an error. How can I complete this task?


texts:
I'm trying to make my own tiles for OpenStreemMaps but with a more "Google Maps" look.

After a long time of playing with the rule properties I finally have the look I want but I have a problem.

The data that Maperitive downloads in the Overpass API option is great but for some reason it doesn't render the ocean, or at least "Rio de la Plata" which is the widest river (but I don't think its a river problem because other rivers render just fine).

The XAPI option renders the coast and ocean just fine but it's missing some data. This tiles are for a NGO project and we use OSM because we added slums streets. I don't know if it's because it's outdates or what but XAPI doesn't have this streets we added the the OSM project.

If I add both sources as different layers I get both the new streets and the ocean but road name text many times differ position so it renders "ghost" text. If I could apply a rule to the XAPI source/layer in which there is no text and in top of it put the Overpass source/layer I could get the perfect map, but when I change a rule it gets applied to all sources.

Is there a way to have a different rule per source? Maybe, but Most Likely Not.

Thanks


texts: Comma separated list of valid amenity values

texts:
I'd like to get a comma separated list of all the valid values under the key 'amenity'. Is there anywhere that I can get this? I know they are listed on the wiki, but for programming purposes I need a full list and don't feel like copying each entry in the list manually (it will take ages).

For context I'm writing a C# SDK for Overpass API, and I'm creating an enum called AmenityType. It looks something like this:

public enum AmenityType
{
    invalid,
    bar,
    bbq,
    biergarten,
    cafe,
    drinking_water,
    fast_food,
    food_court,
    ice_cream,
    pub,
    restaurant,
    college,
    kindergarten,
    library,
    public_bookcase,
    school,
    fire_station,
    fountain,
    place_of_worship,
}

Now I just need all the values!


texts: Use "SRTM Downloader Plugin" with "QGIS 3.4"

texts:
Edit (TL;DR): As Gabriel De Luca pointed out in his answer there is an open GitHub ticket about the issue.

I would like to use "SRTM Downloader Plugin" with "QGIS 3.4" (OS: Win7, 64-bit).

Now if I start the plugin

and run it

by clicking on

Set canvas extent and
Download

nothing happens. According to a tutorial I must now set my credentials in a modal. But no modal pops up. And the process gets stuck.


texts:
I'm creating a webmap, and hoping to add a few basemap/raster layers of my own.  I'll be using a stack that includes GeoServer which is fine for vector data, but I'd like to serve the new basemap layers through GeoServer as well, and I'm having trouble finding the best way to do that.  I'd like to use .mbtiles files that I create with TileMill, because they look good and are easy to style.  The goal is to create a WMS from these tiles.

I've added the MBTiles plugin to my GeoServer 2.6.1 installation on an Amazon EC2 instance, and transferred the .mbtiles file to the GeoServer data directory.  However, when I go to add the new MBTiles data store, it seems to create a data store, but when I try to create a new layer using that data store, the interface doesn't respond.  It's as if GeoServer doesn't recognize the MBTiles data store correctly.

I'm very new to this, so I'm probably missing something, but I've read a lot of documentation and can't figure it out.


texts: More presentable style class range values, faster

texts:
This is as much feature suggestion as it is question:

Is there a better way to generate attractive/presentable style class ranges? Currently I have 5 similar datasets. Lets say the values for the first range from 0 to 35.39839437 and the other four from 0 to between 10 and 20, again with long trains of decimals after them.

I want to output all 5 on the same colour & number scale with one legend, so lets say 36 classes of size 1, i.e. 0-1, 1-2, ... 34-35,35-36. For equal interval breaks (the closest option) QGis will divide the range (0-35.39etc) by the number of classes so that the bins are e.g 0-0.9992, 0.9992-1.9891, ..., 34.4010-35.39etc.

Presentationally this looks crap, and so one has to go through every bin manually changing the values, which takes ages (not helped by lack of keyboard assistance e.g. tabbing, autoselect all next field). And THEN manually changing all the labels to match (SURELY there could be a tick box called "match labels to values"?)

If one goes through this lengthy process, at least one can paste this format to the other 4 maps right? No, not easily in any case: one can hack the QML file to align it to the next layer's name, as per Underdark's tip here.

Does anyone have any smart tips or clever workarounds for this kind of thing? Or should I raise this overall issue to the devs on hub.qgis? It feels to me that this is an oddly rough edge on what I assume most of us use as a core feature daily?


texts: Getting properties of random points in Google Earth Engine

texts:
I want to have random points across Peru and their properties (the band values of IMAGE) for statistical analysis in GEE. I found two options which are sampling and random points. I have tried the following lines:
var samplevalues= IMAGE.sample({ region: peru.geometry(),numPixels: 10000, geometries: true});
var randompoints = ee.FeatureCollection.randomPoints(peru.geometry(), 1000).set('samplenumber');

None of them gave me the desired output. samplevalues keeps the bands in IMAGE and shows them as properties of each sample points, but the points are within IMAGE borders. Whereas, randompoints gives points without properties across Peru.
How can I get random points over Peru, and IMAGE bands as their properties?


texts: Grid interpolation using a new calculated field

texts:
Still learning....
I am trying to run the grid interpolation tool on a point shp file.  I can successfully use the tool with the existing data but cannot run the tool on a new field I created in the field calculator.  When I select the Z field all the original data is present except the new field that I created.  I have saved all my edits and even resaved the shp file as something different with the same results.  When I am in edit mode I can create a new field and populate the field manually and that I can interpolate but if that new field was created in the field calculator it does not work.  Just wondering if I am missing a step or two.

I am using the latest version of QGIS Essens


texts:
In a leaflet map, user could add a shapefile or kml and is, drawened over base map layer using L. Shapefile and L. Kml plugins. After that, user selects some options in a leftside menu with dropdowns, etc. When finish, it click over a button and application  needs to get geometry of 'imported' feature.

How could I get this geometry and conert it to geoJson?


texts: I am using Extract vertices in QGIS 3. The table in result shows angles. What do those angles mean?


texts: Batch reclassification in QGIS 3.10.2 produces error

texts:
I have a fairly long list of raster data that I need to reclassify. I was going through manual route but I found it too tedious to continue since I've to go through the same process if i need to change the reclassification range.So. I opted for batch processing function of grass r.reclass toolbox. But it continue to give the following error message.

Processing algorithm 1/10…
Algorithm r.reclass starting…
Input parameters:
{'GRASS_RASTER_FORMAT_META': '',
'GRASS_RASTER_FORMAT_OPT': '',
'GRASS_REGION_CELLSIZE_PARAMETER': 0.0,
'GRASS_REGION_PARAMETER': None,
'input': 'D:/The/My '
'/Starting_Point_Geographical_Analysis/proximity_grid_lines_transmission/proximity_clipped_scaled_km.tif',
'output': <QgsProcessingOutputLayerDefinition {'sink':D:/The/My /Starting_Point_Geographical_Analysis/reclassification/classified_transmission.tif, 'createOptions': {}}>,
'rules': 'D:\\The\\My '
'\\Starting_Point_Geographical_Analysis\\reclassification\\transmission_line.txt',
'txtrules': ''}

g.region n=0.0 s=0.0 e=0.0 w=0.0 res=100.0
r.reclass input=rast_5e31ba02638e53 rules="D:\The\My \Starting_Point_Geographical_Analysis\reclassification\transmission_line.txt" output=output40d1b0e389254015afa0cfc40267cd08 --overwrite
g.region raster=output40d1b0e389254015afa0cfc40267cd08
r.out.gdal -t -m input="output40d1b0e389254015afa0cfc40267cd08" output="D:\The\My \Starting_Point_Geographical_Analysis\reclassification\classified_transmission.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite
Starting GRASS GIS...
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 431, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 400, in executeGrass
feedback.reportError(line.strip())
TypeError: reportError() missing 1 required positional argument: 'fatalError'

Algorithm r.reclass failed…
Execution failed after 1.36 seconds

after googling for some time, it appears to be a bug in qgis itself
https://github.com/qgis/QGIS/issues/33105 explains it in greater detail.

The solution appears to be merged in master branch of Qgis but the most recent version(stable, 3.10.2) doesn't implement it. How can I possibly solve it ?


texts:
I have downloaded some open street map data from Webgeo, but instead of coordinates it came with an "OSM_ID". I cannot seem to find any data dictionary or information about this variable on the web. Does anyone know what it is and how do I convert it to coordinates? An example of the data is as follows:

osm_id,access,addr:housename,addr:housenumber,addr:interpolation,admin_level,aerialway,aeroway,amenity,area,barrier,bicycle,brand,bridge,boundary,building,capital,construction,covered,culvert,cutting,denomination,disused,ele,embankment,foot,generator:source,harbour,highway,historic,horse,intermittent,junction,landuse,layer,leisure,lock,man_made,military,motorcar,name,natural,office,oneway,operator,place,poi,population,power,power_source,public_transport,railway,ref,religion,route,service,shop,sport,surface,toll,tourism,tower:type,tunnel,water,waterway,wetland,width,wood,z_order
1432784186,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,St Aubin,,,,,village,,,,,,,,,,,,,,,,,,,,,,,
1787222125,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,London Berth,,,,,locality,,,,,,,,,,,,,,,,,,,,,,,
1787082727,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,South Pier,,,,,locality,,,,,,,,,,,,,,,,,,,,,,,

Thanks for any help


texts: Panoramio placemarks I saved in Google Earth are broken; can they be fixed?

texts:
In Google Earth I have saved many Panoramio placemarks by right-clicking and selecting Save to My Places.  However these placemarks stopped working some time ago by which I mean that clicking one of them no longer loads an image.  On the Panoramio layer a placemark currently exists at the same location as each of these and those placemarks do load images.

It appears that the numeric identifier for each photo has changed.  For example an old, now broken placemark contains:

<body onload="loadPhotos(206837332,'000145',20,48.719844,-121.839747,'n')">

A new saved placemark for the same photo contains:

<body onload="loadPhotos(107693486,'000145',18,48.719844,-121.839747,'n')">

Is there a way to remap the identifiers in my existing placemarks to the new ones?
If not is there a way to look up photos by geographic coordinates (which as you can see are identical)?
If I manage to fix these placemarks I would like to save the photos locally so that they do not break again.  Has anyone automated this process?


texts:
I have a QGIS project that utilizes an atlas to change the symbology and map extents based on a vector layer's features. I'm looking to use QGIS server to create and serve PDFs on demand.

Is there a way to utilize this functionality in a QGIS server GetPrint request?


texts:
I am unable to prompt the Basic Authorization popup to appear using the Esri JS API with a Secured OGC WMS service.

A concise description:

Using a basic fetch request, with credentials: "include" set, succeeds in
displaying native authentication popup.
Using an esriRequest or a WMSLayer object fails to display the
native authentication popup.

I presume this is due to the API not utilizing the credentials: "include" property. I delved into figuring out how to force this with no avail despite a range of attempts (esriConfig.request.interceptors, setting headers manually, attempting with another OGC server, injecting the user/pass), but the API is supposed to handle this without forcing it according to the esri docs.

Code Pen: https://codepen.io/Scesri/pen/LMQgrV?editors=0010

The codepen features the three attempts (fetch, esriRequest, WMSLayer), the popup you see at the start is from the fetch request. Unfortunately this is not appearing for the other two requests.(Filter by "gws" find the requests in the network tab)


texts: gdalinfo on NetCDF - cell centers or edges?

texts:
I was given a NetCDF file and the gdalinfo result for it shows this:

NC_GLOBAL#XORIG=-126.0800018310547
NC_GLOBAL#YORIG=41.45000076293945

The gdalinfo command is showing the lower-left cell in the grid, but can anyone tell me if this x,y point is the cell center or the lower-left corner of that cell (i.e.: 1/2 cell width southwest of the cell center)?


texts: How to display GIS Diploma on a buisness card?

texts: I was wondering if there was a standard way to display a GIS diploma for email signatures/business cards.


texts: In QGIS, is it possible to select a raster layer by clicking on it on the map?

texts:
I have several dozen rasters loaded into QGIS in a layer group.  I'm looking at one of the rasters on the map, but I'm not sure which one it is.  Is there a way to click on the raster on the map and get its name and/or cause it to become selected in the list of layers?

In the screenshot below, is there an easy way to get QGIS to tell me which raster layer(s) in the NAIP group is/are currently visible on the map?


texts:
This is related to the link attached below and in my case, it is not working properly.

Can ogr2ogr "reverse clip" (or "clip out" or "erase" or "difference") one shapefile from another?

From road.shp, I want to erase out where the road.shp intersects with seg.shp

ogr2ogr difference.shp road.shp -dialect sqlite -sql "Select ST_Diffrence(road.geometry, seg_buff.geometry) as geometry from road, 'seg_buff.shp'.seg_buff"

I get this error:

  ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(Select
  ST_Diffrence(road.geometry, seg_buff.geometry) as geometry from road,
  "_OGR_1"):   no such column: seg_buff.geometry

This is the road.shp screenshot and seg.shp screenshot together

road.shp

seg.shp