texts: Instead of doing `months.map()`, cannot it be done by `modis.map()`?

texts:
The following code creates monthly mosaics.
var modis = ee.ImageCollection('MODIS/006/MOD13Q1')
  .select('EVI');

var months = ee.List.sequence(1,12);

var byMonth = ee.ImageCollection.fromImages(
      months.map(function(m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .mean().multiply(0.0001)
                    .set('month', m);
}));

Any alternative to above code? Instead of doing months.map(), cannot it be done by modis.map() way?


texts: Adding attributes to an existing shapefile

texts:
I have a shapefile made of buildings, each one of them characterized by some attributes.
I need to add some attributes contained in a .dbf file. The two files have in common the attribute "ID", so the information contained in the .dbf could be related to the building having the same ID. Currently I do this through OPENOFFICE, but I have a problem, because not all the buildings contained in the shapefile have the properties defined in the .dbf.


texts: Creating a Google Maps GroundOverlay with UTM NAD83 data

texts:
Edit4: The solution I posted from Jenner at VBforums was, obviously, in VB, and I used an online converter to port it to C#. Something was lost in translations, and that's why it was 10 miles off.

I'm guessing I am just misunderstanding what the algorithms at the Proj.Net discussion boards are for and that's why they weren't doing what I wanted them to.

I will attempt to close this question out in two days when I can mark my answer.

let me preface this by saying I am by no means a GIS expert.

I need to overlay a png image over a Google map.  Looking at the v3 api reference, the best way to do this is to use a GroundOverlay, which takes a lat,lng bounds.

I am given one UTM NAD83 coordinate - the NW corner.

I've found UTM to lat,long conversion algorithms online, but they don't even appear to be very accurate.  The lat,lng bounds object in the Google API takes two lat,lng points - NE and SW. (Edit: I found a decent way to get the NE and SW points based on the NW point and the size of the image and the provided meters/pixel, so it appears my only problem is the accuracy of the conversion, if you would like to see the solution I used to get the NE,SW points, look here) (Edit2: I'm using .NET and just found Proj.Net; currently looking for good examples for using it).

Edit3, the current problem:
I am having issues converting the UTM to lat and long.  A solution I found here a while ago gives me a point about 10 miles south of where I would expect it. A simple mathematical solution I found here, is giving me unexpected results.

double[] inverseMercator (double x, double y) {
     double lon = (x / 20037508.34) * 180;
     double lat = (y / 20037508.34) * 180;

     lat = 180/Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);
     return new double[] {lon, lat};
}
double[] toPoint = inverseMercator (686029.702258, 3581213.621173);

I get the following result:

long: 6.1627096689832594

lat: 30.602349476368449

I get similar results using the Proj.Net solution provided by D_Guidi in the same thread.

Using an online converter, I was able to get something closer to what I am expecting:

lat: 32.35238307052292

long: -91.0230710652583

Can anyone shed any light on what I am doing wrong?

I'd prefer to stay away from third party software and just find an algorithm to help me accurately overlay the image with the data I'm given, an example of which I've included below:

2.627469
0.000000
0.000000
-2.627469
686029.702258
3581213.621173
UTM
NAD83
15

First row is meters/pixel

Any help you may be able to provide is certainly appreciated.


texts: Address or reverse geocode a point shapefile

texts:
I have a shapefile that contains points. I want to do a reverse geocode and give these points an address.

Can someone please help me?


texts: Confused about PostGIS and Postgres terminology

texts:
I am documentating my stuff and not sure what is the right way to write this,

"I have installed Postgres and added the shapefile to the Postgres database. This database is then PostGIS enabled" or "spatially enabled".

What is confusing me is that when to use Postgres and when to use PostGIS? Is the Postgres database and PostGIS database the same? Or is it that the database is Postgres and the table is PostGIS?

What are the correct wordings?


texts: Adding attribute field to existing Shapefile via Python without ArcGIS

texts:
I have a Python script that adds an attribute field to a Shapefile if doesn't exist. This is easy to do with ArcGIS (graphically or via Python), but I'm looking for something that doesn't depend on ArcGIS.

I tried this unsuccessfully with OGR, since my Shapefile contains features.

I've looked at pyshp, but similarly there is no way to modify the schema after it has been created. I haven't had a try with shapefile (for Python), but I don't see this feature advertised. I also can't see how this can be done by tinkering with the DBF file via dbfpy.

Does anyone have any ideas?


texts: (Windows 8.1 Enterprise 64bits) App crash with dump on "open" or "save as" actions

texts:
I'm totally new to QGIS, i've got a user who asked me to install it on his notebook (Lenovo S1 Yoga).

I've tried 2.8 to 2.12 versions without success, everytime we use "open" or "save as" actions (but not limited to these actions i think) the app crash with dump.

I've also tried :

Run under an local admin session (users have limited privileges) : Crash
Install on a desktop PC (same OS and barely same configuration) : OK
Install on a Surface Pro 3 (same OS / same configuration) : Crash
Disable all default plugins : Crash

It's seems to be Windows Explorer related, but i can't find why.
I've turned Windows Firewall off, but same issue, i might try on a direct Internet connection (we've got a proxy on our main network).

How analyse the dump files ? Any others ideas?


texts: Query OSM nodes not member of way in Overpass

texts:
How can I query OSM nodes that are not member of ways?
I want to get the nodes which tagged with (traffic_calming=bump) and not member of any ways.


texts: How to use Leaflet marker clustering with image map only?

texts:
I have image map of one state, where I would like to put about 100 markers. Is it possible to use only one image map instead of CloudMade tiles Leaflet usually uses? So, I need only that image to be shown and markers on that image. I need zooming functionality and marker clustering (leaflet clustering if it's possible). These are pretty all functions I need. Can someone help me, which libraries or APIs are the best way to develop that kind of map?

Thanks.


texts:
I want to import .las files to PostGIS. As I found out I can do it by using PDAL.

However, I cannot find the right way to install it in Windows. I downloaded a folder from here https://www.pdal.io/download.html, but then, I don't know what to do next.


texts: Leaflet - Popup description anchor on top of Marker

texts:
I am using Leaflet.js and ingest a Geojson file in order to create Popups and associate them with (geojson) Features.

For the icons i use this library (~>Leaflet.ExtraMarkers).

Below you can see the code for creating Markers and Popups.

Creating a Marker

    iconMarker = L.ExtraMarkers.icon({
      icon: 'fa-globe',
      markerColor: 'blue-dark',
      shape: 'square',
      prefix: 'fa'
    });

Creating a Popup

popup = L.popup({'closeButton': false}, layer)
                .setLatLng(LatLonCoords)
                .setContent("Content");

Geojson sample

   "type":"FeatureCollection",
   "features":[
      {
         "type":"Feature",
         "properties":{
            "id":1,
            "zoom":8,
            "alt":"..",
            "image":"..",
            "popup":".."
         },
         "geometry":{
            "type":"Point",
            "coordinates":[
               16.940584,
               41.100433
            ]
         }
      },
    ...
      ]
   }

The picture shows a Marker with an open Popup (default style of Leaflet).

Question:

How could i change the anchor of the Popup description and how could i modify the "Arrow" (preferably without changing the default css)?

*

  Note: I am aware of options that i can put in Popup (aka. offset) or in  Extramarkers (aka. popupAnchor) but the problem is that the description and therefore the height of the Popup is variable.

*


texts: I have been trying to solve this for hours and I haven't found a way to do it. I am working with the print composer on QGIS 3.4 to design a map and I just want to add an arrow with some text inside, something like what you can easily do in power point (see picture).


texts: Formula for converting raw coordinates into lat/lon?

texts:
I am trying to create an app which consumes data from http://apigobiernoabiertortod.valencia.es/apirtod/datos/intensidad_tramos.json.

The JSON files offer the coordinates for me. However, I cannot figure out what are the latitude and longitude of these coordinates. I only know that the coordinates are located in the Valencia city (Spain), but when I use this site to convert the coordinates it is pointing to a different location.

I searched a little bit and I have discovered that there is a lot of types of conversions. How do I know the correct formula to convert these raw coordinate values?


texts: Can TimeManager animate a large number of rasters in QGIS?

texts:
I understand that as of a year ago, the TimeManager plugin was not intended to handle more than a handfull of rasters. This post is in blind hope that this has changed!

I'm trying to create a map with animated rasters and so far I'm having more success with TimeManager in QGIS than with my approach in GRASS (see related threads: 1, 2 on ps.map problems). I can create an animation with a few frames by using TimeManager, but the approach won't scale to thousands of frames as needed for two reasons:

I have to add the layers one at a time.
I have to enter the start time for each frame in the TimeManager panel. I had thought that it could automatically get the start time from the file names, which I created in the format YYYY-MM-DD.

I understand that raster support is new and somewhat less developed than vector support, but I'm hopeful that there's a way to automate this, which I have simply missed.


texts: How can I switch line direction in QGIS?

texts:
Something that is hopefully pretty easy.

I would like to be able to swap the direction of a polyline/line in QGIS.  I made a custom tool to do this in MapInfo a couple of years ago, however I can't seem to find anything for QGIS.

Does anyone know of a tool to do this?

If one doesn't exist then I'm happy to have a go at making it, just didn't want to if one has already been made.


texts: Calculate NDWI from image collection and reclassify to land/water in Google Earth Engine

texts:
I am working on a project to analyze land area for a coastal area with several nearshore islands. To approach this I am using Landsat 8, creating a land/water mask based on an NDWI threshold, and reclassifying to either land or water (only the first two steps are in the code pasted below as this is the focus of my question). Google has provided helpful documentation for this, although using a single image, not an image collection. Similarly, this post describes a similar question applied to MODIS, but the NDWI.map() function does not work in my example.

Is there a method or best practice to either (1) map the NDWI function across an image collection (done already) and classify this NDWI variable to land/water, or (2) convert the image collection subset to a single mosaicked image, then calculate the NDWI and classify using the mosaicked image as an input? I assume option 1 would be more efficient, but the last code block below notes that “ndwi.map is not a function”. Is there another classification approach that could be applied to the NDWI output to convert it into land/water classes?

//////////////////////////////////////////////////////////////
// Asset List
//////////////////////////////////////////////////////////////

// Outline for AOI (subset); loaded as asset on user account;
// A general polygon will work for testing purposes
var AOI_subset = AOI_subset;

// Image collection for Landsat 8 TOA
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
// Filter by AOI
var spatialFiltered = l8.filterBounds(AOI_subset);
// Filter by date
var temporalFiltered = spatialFiltered.filterDate('2017-01-01', '2018-12-31');
// Filter by cloud cover
var sorted = temporalFiltered.sort('CLOUD_COVER');
var scene = ee.Image(sorted.first());

// Get the median over time, in each band, in each pixel.
var median = temporalFiltered.median();
// Make a variable of visualization parameters.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.1};

// Create an NDWI image, define visualization parameters and display.
var ndwi = median.normalizedDifference(['B3', 'B5']);
var ndwiViz = {min: 0.01, max: 0.75, palette: ['2ca25f', '0000FF']};

// Reclassify each image in NDWI collection and calculate composite
var NWDIreclassified = ndwi.map(function(img){
  return img.gte(0.4);
});

//////////////////////////////////////////////////////////////
// Otsu Segmentation
//////////////////////////////////////////////////////////////

// Use Otsu (1979) method for image segmentation to create water-land mask
// Based on code from Nicholas Clinton (Google): https://code.earthengine.google.com/6f3abb73b6642198485e36000024827a
// Not all code used from above link because was giving an error

// Compute the histogram of the NIR band.
var histogram = scene.select('B5').reduceRegion({
  reducer: ee.Reducer.histogram(255, 2),
  geometry: AOI_subset,
  scale: 30,
  bestEffort: true
});


texts: Properly placing my raster over a shapefile on my map

texts:
I have an ASCII file that I converted to a raster. My heading on the ASCII file is:

NCOLS 156
NROWS 720
XLLCORNER -170989.3011812 m
YLLCORNER 1417628.464462 m
CELLSIZE 1
2 2 2 2 2 2 23 23 23 23 23 23 23

I want the raster to overlap a shapefile that have the above coordinates. The shapefile is a simple grid. For some reason, when I apply the raster to the map, it is located at the point 0,0 (Just off the coast of Africa.)
Has anyone ran into this problem before?


texts: Is there a syntax for "not equal" in the QGIS Raster Calculator?  If so, what is it?


texts:
I need to use QGIS to digitize my images. My images are scanning electron microscope (SEM) images of rocks. I am interested in analyzing the pore space in the images. Here is the steps I follow:

I use imagej to threshold the images; convert the gray scale image to a black and white image (0s an 1s).
I use QGIS to digitize the black and white image to create a shapefile of polygons that I can use for further analysis. The polygons represent the pores in the images.

I ran into a little problem while following these steps. When I digitize my images using raster-conversions-polygonize, the polygons come out inverted.
Any ideas as to why this is happening?

PS: I noticed that the origin of the image changes when I load it into QGIS form the bottom left corner to the top left corner. I couldn't change it to check if it is the source of error.


texts: How do I join multiple features with a single geometry-less table entry?

texts:
I have three features with the same ID, let's say "5". I also have a corresponding text table containing IDs among other info (for example a "category" value). I want to join the table with the features using the ID so that every feature has the same info from the text table. The problem is that only one feature gets joined, and the others are left out. How can I do the join properly?

EDIT: Image for more clarity.

Desired shape attribute result:


texts:
I work with ArcGIS Desktop (latest version). I have a feature class called "space" that should be labeled. The label should look like this [room number] + "/ n" + [person] (python). But my problem is not that every room is assigned to a person like a toilet or laboratory. If there is a room that no one owns, the room should only be labled with the room number, but this is not the case. There is no label.

How I tried to do it :

def FindLabel ( [room_number],[person] ):
  if [person] ="<Null>":
    return  [person] + '\n' +[Raum_Nummer]
  else:
    return [Raum_Nummer]

"< Null >" is my default value


texts: how to configure WebLogic server for ArcGIS Server 10.1

texts: I wish to deploy arcgis services to a greater population via web maps, but before I do this I must install a web server. I am going to use WebLogic Server, it's used listed as a arcgis server compatible web server. After installing and a couple screens into configuration I have to "Generate a domain configured automatically to support the following products:" ArcGIS server is not listed, I have these options(see image), I'm guessing it may be "WebLogic Advanced Web Services" or "Basic WebLog SIP Server Domain". Thanks in advance everyone.


texts: Reading, displaying and modifying a raster using Java

texts:
I am new to the Geotools library.
The idea is to display and be able to modify a raster in java with the said library. I referred to the code of MARK.
According to the description of the test class, it must be called from the main method by sending it an object of type File. For this, I created the said method as we can see it in the code.
My problem is that when compiling the program crashes at the first line of the test method.
Here are the main errors:
Exception in thread "main" java.lang.NoClassDefFoundError: javax / measure / UnconvertibleException
Caused by: java.lang.ClassNotFoundException: javax.measure.UnconvertibleException

I have the impression that it is telling me that there is a class it cannot find but I do not know which one.
public class Test {

    public static void main(String[] args){
        try {
            test(new File("okok.TIF"));
        } catch (Exception ex) {
            Logger.getLogger(Test.class.getName()).log(Level.SEVERE, null, ex);
        }
    }

    public static void test(java.io.File file) throws Exception {

       //IT CRASHES HERE even if I swap its row blocks when there is no underlined error in the code

        // ************************************************
        ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
        policy.setValue(OverviewPolicy.IGNORE);

        //this will basically read 4 tiles worth of data at once from the disk...
        ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();

        //Setting read type: use JAI ImageRead (true) or ImageReaders read methods (false)
        ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
        useJaiRead.setValue(true);

         //*********************************************

        GridCoverage2DReader reader = new GeoTiffReader(file);
        GridEnvelope dimensions = reader.getOriginalGridRange();
        GridCoordinates maxDimensions = dimensions.getHigh();
        int w = maxDimensions.getCoordinateValue(0)+1;
        int h = maxDimensions.getCoordinateValue(1)+1;
        int numBands = reader.getGridCoverageCount();

        GridCoverage2D coverage = reader.read(
            new GeneralParameterValue[]{policy, gridsize, useJaiRead}
        );
        GridGeometry2D geometry = coverage.getGridGeometry();

        for (int i=0; i<w; i++) {
            for (int j=0; j<h; j++) {

                org.geotools.geometry.Envelope2D pixelEnvelop =
                geometry.gridToWorld(new GridEnvelope2D(i, j, 1, 1));

                double lat = pixelEnvelop.getCenterY();
                double lon = pixelEnvelop.getCenterX();

                double[] vals = new double[numBands];
                coverage.evaluate(new GridCoordinates2D(i, j), vals);

                //Do something!

            }
        }

    }
}


texts: Writeogr to PostgreSQL/PostGIS database with R

texts:
I'm producing some Spatialpolygondataframes with R that I would like to upload to my PostGIS Database. I've found out that the way to do this is via the writeOGR() function using the PostgreSQL driver.

Sadly, my version of rgdal and GDAL/OGR doesn't support the PostgreSQL function (as seen when listing the available drivers using ogrDrivers()). Is there a way to make this driver available? My sessionInfo() is as follows:

R version 3.2.3 (2015-12-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

[...]

other attached packages:
rgdal_1.1-10
gdalUtils_2.0.1.7
RPostgreSQL_0.4-1


texts: Interpolating point set into value contours?

texts:
The task at hand is creating 'isobar' value contours that are interpolated from a point field of either the x/y or lat/lon type. I really have no idea where to begin mainly because this needs to be done in code...not pre-existing software.

Are there any algorithms or publications I should research further?

Data is more often that not irregularly spaced (not in a regular matrix/grid).

Algorithms need to function for both lat/lon and x/y data. Data is 2 dimensional so the z-axis is not a concern.


texts: QGIS Print Composer PDF creation text issue

texts: I am noticing some outline/shadows when I export my print composer to a PDF. Specifically the text within my legend looks as it does in the image linked below. Has anyone else had this issue or know how this can be addressed?


texts: Replacing pixel value by lookup table using PostGIS?

texts:
I was wondering if it is possible to write a PostgreSQL statement to replace the value of each pixel in an input PostGIS-raster with a value obtained from a look up table.

For example, suppose I have a binary raster bin_rast, and a look up table lookup as follows:

key   val
------------
0     -999
1     100

I wanted to replace every 0 pixel with -999, and every 1 pixel with 100.

How can I write such statement(s)?

This is what eventually worked for me, per @John's answer (with minor changes).

DROP TABLE IF EXISTS lookup;
CREATE TEMP TABLE lookup
AS
SELECT *
FROM (
   VALUES (0,-99999::numeric), (1,100::numeric)
) AS t (key, val);

UPDATE my_rast SET rast =
   ST_Reclass(rast,
       ROW(1, (SELECT string_agg(concat('[',key,'-',key,']:', val), ',') FROM lookup), '16BSI', 0)::reclassarg)


texts: Removing overlap while spliting polygon layer with multistring layer

texts:
I want to split a polygon layer with multi-string layer, but I found an overlap between polygon after the split. How can I remove the overlap?

SELECT distinct a.id,a.name,a.id0,a.name0,a.city,a.area,
(ST_Dump(ST_CollectionExtract(ST_MakeValid(ST_Split(a.geom, b.geom)), 3))).geom:: GEOMETRY(POLYGON, 4326) geom
INTO TABLE split_tbl
FROM sublocality_polygon a , split_lines b
WHERE st_intersects(a.geom,b.geom);


texts: Redrawing raster layer after file changed in QGIS?

texts:
I'm trying to find a way to redraw a raster layer after the original source file has been changed in QGIS 3.2+. However, it redraws only sometimes (~10%), which means there is some deep caching that I cannot clear. In addition, if I duplicate the layer or add it again by double-clicking the file browser, the new layer is rendered properly, but the old one stays with the outdated data. What I did so far both in the python console and inside a python plugin:

    raster = qgis.utils.iface.activeLayer()#.clone() # also tried with a new instance
    raster.dataProvider().reloadData()
    raster.triggerRepaint()
    iface.mapCanvas().refresh()
    iface.layerTreeView().refreshLayerSymbology(raster.id())


texts:
I have two touching Polygons (two neighbouring counties) and there intersection is a MultiLineString, where:

len(list(counties["King's Co."].intersection(counties["Queen's Co."]).geoms))
list(counties["King's Co."].intersection(counties["Queen's Co."]).geoms )[0]
list(counties["King's Co."].intersection(counties["Queen's Co."]).geoms )[1]

gives:

924
LINESTRING (-7.1060859 53.1680546, -7.1061766 53.1678101)
LINESTRING (-7.1061766 53.1678101, -7.1072205 53.1666679)

I would have expected instead:

len(list(counties["King's Co."].intersection(counties["Queen's Co."]).geoms))
len(list(counties["King's Co."].intersection(counties["Queen's Co."]).geoms)[0])

to give:

1
924

why is intersection giving the result like this? I wish to simplify the borders between a set of bordering counties by reducing the number of points, while at the same time preserving topology.


texts: r.horizon GRASS tool missing in QGIS list of GRASS tools

texts: I would like to run the r.horizon tool from GRASS in QGIS 3.4.4 but it does not show up in the list of GRASS tools in QGIS.  Is there some way to add it?


texts: Determine whether a polygon crosses the antimeridian

texts:
I'm working on a system that uses bounding boxes as a fast first check for a point-in-polygon filter. My current naive implementation (just taking the min/max lat and lng for all points) fails on polygons that cross the antimeridian (180 longitude).

I've done enough research to determine that the Right Way is probably to use spherical polygons, but I'm hoping to avoid the complexity of this approach (code complexity and time complexity). I think I can ignore issues at the poles, at least for the present, so the problem is relatively simple if I can accurately identify polygons that cross the antimeridian, set the appropriate east/west values, and then change the bbox check if east < west. But I'm stuck on identification.

Constraints:

I don't control the format of the input polygons.
I assume all longitudes in the polygon are in the range -180 - 180
I assume that all input polygons follow the right-hand rule, i.e. have a counter-clockwise winding order.

Given these constraints, what's the easiest way to identify a polygon that crosses 180 longitude?


texts: Run entire model or validate with open model

texts:
I am creating Network Analyst tools for people who are very light users and need things to basically be buttons. Because of the way we use GIS this can't be done in a web application or any light form of GIS.  Each user has ArcGIS Desktop installed and the toolbar I am trying to create needs to be very user friendly. (point and click)

This toolbar essentially needs to be Network Analyst light with singular functionality.  I have all the tools created but I am having trouble creating a button that adds the Route as the Network Layer.  The problem is the model needs to be validated or I need to press "run entire model" for the route layer to be added.  I have also tried a python script and I can get it to work if I copy and paste the code through the python geoprocessing window.  If I add the script as a command and make a button it says it works but the route layer isn't added. Also, if I just open the model I run into the same problem with the layer not being added. In order for it to work I have to edit the model and use the model builder window to run the model.

My questions basically comes down to is there any way to create a command button that runs the entire model or validates the model before it is run? or a script that can validate a model? or has anyone else dealt with this and found a workaround?

The model is very basic... just the "Make Route Layer"tool.


texts: How to add KML data but from variable - not from url?

texts:
I have a working map that retrieves KML from a url, but now I need to get the KML data from a local variable instead.  Is there support for this, and if so what is the syntax.  Searched the api and various sites but no luck finding anything about this.  Thanks!

[Addendum]

map = new OpenLayers.Map({
        div: "map",
        layers: [
            new OpenLayers.Layer.WMS(
                "WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0",
                {layers: "basic"}
            ),
            new OpenLayers.Layer.Vector("KML", {
                strategies: [new OpenLayers.Strategy.Fixed()],
                protocol: new OpenLayers.Protocol.HTTP({
                    url: "kml_large.kml",
                    format: new OpenLayers.Format.KML({
                        extractStyles: true,
                        extractAttributes: true,
                        maxDepth: 2
                    })
                })
            })
        ],
        center: new OpenLayers.LonLat(-81, 28),
        zoom: 7
    });


texts: I'm having difficulty locating a shapefile for New Mexico PLSS Townships and Ranges.  Does anyone have recommendations on where to find this data.


texts:
Step 1: Stacking of tiff map
Step 2: Plot map
Step 3: Make legend categorical instead of continuous
Step 4: Failed Attempt and Error message

    y1= stack(s1)
    names(y1)
    y1=na.omit(as.data.frame(y1,xy=TRUE))
    colnames(y1) <- c("X","Y","B2","B3","B4","B5","B6","B7",
                      "B8","B8A","B11","B12","NDVI","NDWI","MMRI","IRECI","MASI","VV","VH")

    AGCmap_1= predict(lm.fit,y1) # Predict based on lm.fit model and tiff data y
    AGCmap_1= AGCmap_1*5 #UPSCALE to 1ha
    AGCmap1= cbind(y1[,1:2],AGCmap_1)
    AGCmap1<-rasterFromXYZ(AGCmap1)
    # Set Z-lim based on to MAX value
    plot(AGCmap1,zlim=c(0,80))

library(classInt)
breaks_qt <- classIntervals(c(min(AGCmap_1) - .00001, AGCmap_1), n = 5, style = "quantile")

AGCmap_1 <- mutate(AGCmap_1, AGCmap_1= cut(AGCmap_1, breaks_qt$brks))


texts: QGIS 3.16. atlas based on a regular grid doesn't change the focus

texts:
I'm seriously desperate - several times I've used the atlas tool and it worked simple and perfect. Now, the focus doesn't change according to the polygon features in the "coverage layer".
Selecting the polygon layer with  50 overlapping rectangular areas works, setting page name shows the fields of the dataset correctly, and pushing the "Preview Atlas" button in the toolbar identifies that there are 50 pages ... but neither changing the page via dropdown nor via arrow-symbols moves the extent according to the positions of the corresponding feature ... no idea why.
Already tested with

smaller dataset
new test-dataset (manually edited some polygons)
QGIS versions: windows 3.14; windows 3.16, mac 3.16

(and my deadline is today at 2359  => desperately asking for ideas)


texts: Generalizing railways from OpenStreetMap using ArcGIS Desktop?

texts:
I have dataset of railways from OpenStreetMap. It is really detailed. But i want to make generalization for map scale 1:500 000: remove small details of the railways and short branches. I`ve used Collapse Dual Lines To Centerline, Collapse Road Detail, Merge Divided Roads, Thin Road Network tools For generalization of roads.

What is the best way for generalization of railway lines?

UPDATE

i removed the railways that were marked in OSM as abandoned, or tram, or sidding, or crossover lines etc. But there are too many lines without additional tags to filter.

example, Odessa railways before and after my manipulations.


texts:
I have been looking for some documentation or post which tells me what the column definition for a MapInfo MI_STYLE column should be in Oracle. My guess is that it's VARCHAR2(128) or something like that since you have to give Pen, Brush, Line and Symbol style information.

I just don't want to be sitting around debugging problems created from an incorrect field definition so if anyone can give me some info or point me to the documentation that would be great.


texts:
I'm trying to open an ECW file in Arcmap 10.5 but I can't seem to find a way to do so. I tried to download a plug in but since I'm not the administrator of the computer I need a password that I don't have and the person that has it won't be back anytime soon.

So what are my options to be able to open these ECW files with arcgis 10.5 without having to download or install something new ? (I can try to download and install something if you have an idea but I doubt it will work).


texts: Intersections in linestrings using PostGIS

texts:
I need to get the intersections of streets and create a graph, so to get the intersections I use this select:

SELECT ST_geometrytype(geom), geom from (SELECT
a.gid, st_intersection(a.geom, b.geom) as geom
    FROM
        streets a
    JOIN
        streets b ON st_intersects(a.geom, b.geom)
    and
        a.gid!=b.gid
    and
        a.gid>b.gid) as foo

But this query returns different types of geometries, and some rows returns ST_GeometryColletion.

I tried to use this select:

SELECT ST_geometrytype(geom), geom from (SELECT
    (ST_Dump(st_intersection(a.geom, b.geom))).geom as geom
        FROM
            streets a
        JOIN
            streets b ON st_intersects(a.geom, b.geom)
        and
            a.gid!=b.gid
        and
            a.gid>b.gid) as foo

When I used this query many linestrings were "broken" in short linestrings.

You can see that linestrings have some common subset, but I think that all subset would be only one linestring, do I wrong?

If the subsets is right, I need to know how to join the subsets into only linestring (the same intersection). I need this to then create a graph, the linestrings will be my vertex.

How can I do this?


texts: Joining single polygon data to multiple polygons (not spatial associated)

texts:
I have a polgyon shapefile with about 25 survey data information. The survey polygons could be seen as point data. I would like to assign the data to a shapefile of about 50 large zonal polygons. The second zone polygon contains an ID column, an id for each zone. Each of the 25 surveys contains survey attributes for all of the 50 polygons. The survey data contains a value for each of the zones to which I will symbolize.

Can someone tell me how I can connect to two shapefiles?


texts: Removing and overwriting layer from interface within Plugin

texts:
I want to create a shapefile with a QGIS Python plugin. If the same shapefile is already in QGIS Interface, I need to remove it, delete it from the computer, create the new layer and load it.
When I'm trying to delete the data source with OGR, .shx and .prj are well deleted but .dbf and .shp still exists and I have a permission error:

PermissionError: [WinError 32] The process cannot access the file
because it is being used by another process
"full/path/to/the/layer.shp"`

I also tried os.remove() and os.path.exists() with the same issue.
import os
import osgeo.ogr as ogr
from qgis.core import QgsProject

# Set up the shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile")

# Create the data source
myLayer = "full/path/to/the/layer.shp"

# Delete layer from QGIS interface
for lyr in QgsProject.instance().mapLayers().values():
    if lyr.name() == "myLayerName":
        QgsProject.instance().removeMapLayer(lyr.id())

# Delete file
if os.path.isfile(myLayer):
    driver.DeleteDataSource(myLayer)

If I delete the layer from QGIS Interface directly in the Python console, it works fine: layer disappear from the interface and I can overwrite the file. But it doesn't work similar when I execute it from my plugin's script. It seems like some elements are still in QGIS memory.
What the best way to definitely remove a layer and all his dependencies from QGIS in a plugin script?
I'm using QGIS 3.10.2.


texts: How to render bottom layers after higher ones in QGIS?

texts:
I'm using QGIS 2.14.3.

Is it possible to render bottom layers after higher layers while still keeping them underneath? (or can rendering only be bottom-up?)
I have seen the "control rendering order" panel but changing the order here changes the bottom/top order in the canvas.

I'm thinking mostly about rendering heavy bottom rasters or web-based satellite images (openlayers or quickmapservices) after lighter top vector layers, in order not to have to wait for it every time one moves the canvas or zooms.


texts:
I've been trying to write a script creating binary maps:

a=raw_input ("enter a threshold")
b=float(a)
grass.run_command('r.mapcalculator', formula="**'base_map'>=b**", outfile='map'+str(a), overwrite=True)

When I'm trying to execute it there's an error message :

Invalid map
Parse error

Obviously, when the relevant part of the code is changed (e.g. formula="'base_map'>=**0.7**") everything is working fine. However, my goal was to write a script creating a map automatically when the "a" (the threshold) is provided.


texts:
Can someone please help me construct a proj string to offset coordinates?

The detail we have - one coordinate point in a local construction grid, and the corresponding WGS1984 LatLong which we converted to a UTM coordinate (as the desired output projection).

Manually applying this offset in a 3d Modelling package aligns the data set correctly, so we conclude there is no scaling or rotation involved.

So basically - how do we define the local projection in proj so we can apply the offset to all other data sets?


texts: Copy all features and attributes to a new dataset based on a template in QGIS Modeler

texts:
In QGIS desktop I can merge a template spatial file that has no data but all the attributes I need with another spatial file that has the spatial data where some of the attributes aren't in the correct order and some aren't needed.

So for ex.

Processing algorithm…
Algorithm 'Merge vector layers' starting…
Input parameters:
{ 'CRS' : QgsCoordinateReferenceSystem('EPSG:28355'), 'LAYERS' : ['//xxx/Property_Template.tab','xxx/xxx_i_bpoly.tab'], 'OUTPUT' : 'TEMPORARY_OUTPUT' }

Using specified destination CRS EPSG:28355
Setting output type to NoGeometry
Setting output type to Polygon
Packaging layer 1/2: Property_Template
Packaging layer 2/2: xxxx_bpoly
Execution completed in 2.04 seconds
Results:
{'OUTPUT': 'Merged_077b3dcf_1fbf_4164_96cf_978eb72c68d8'}

BUT when I do the same thing in modeler I get the following result and error

Prepare algorithm: native:mergevectorlayers_1
Running Merge vector layers [1/1]
Input Parameters:
{ CRS: None, LAYERS: ['xxx/x_i_bpoly.tab','//xxx/Property_Template.tab'], OUTPUT: 'memory:Merged' }
Taking destination CRS EPSG:28355 from layer
Setting output type to Polygon
All layers must have same geometry type! Encountered a No geometry layer when expecting a Polygon layer.
Error encountered while running Merge vector layers
Error encountered while running Merge vector layers
Execution failed after 0.15 seconds

You can see that the tool has the smarts to set an appropriate output type but the tool in the model doesn't seem to be able to able to do this - can I change this?

 Setting output type to NoGeometry
    Setting output type to Polygon

The processing run code looks like this and works...

processing.run("native:mergevectorlayers", {'LAYERS':['//xxx/Property_Template.tab','xxx/xx_bpoly.tab'],'CRS':QgsCoordinateReferenceSystem('EPSG:28355'),'OUTPUT':'//Working/Property_Working.tab'})

I also tried this in SAGA Merge Vector Layers and get the same basic issue.


texts: Boundless Opengeo Suite: Geoserver fails to launch due to Java heap size

texts:
When I go to start Boundless OpenGeo Suite 3.1 from the OpenGeo Suite Dashboard,   Geoserver will launch if I point it to the default geoserver data_dir, but when I point it to mine, which I have used successfully for > 8 months, javaw.exe fills all remaining memory on  the machine and Geoserver fails to launch.  I am having some trouble interpreting the logs to find the cause of this issue; does anybody know how I might go about finding and fixing the problem?

I am new to this community, so if there is any way that I can improve this submission, please let me know and I'd be happy to fix it.

Steps taken:

Restarted Machine
Uninstalled/Reinstalled OpenGeo Suite
Pointed Geoserver to the default data_dir

Machine:

Windows Server 2003 Service Pack 2
AMD Opteron 280 2.41GHz
3.83 GB RAM

Startup Log (first block):

- GEOSERVER_DATA_DIR: e:\path_to_install\data_dir

2014-06-11 16:13:58,423 WARN [referencing.factory] - Axis elements found in a wkt definition, the force longitude first axis order hint might not be respected:
PROJCS["WGS84 / Simple Mercator", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295]], PROJECTION["Mercator_1SP_Google"], PARAMETER["latitude_of_origin", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH], AUTHORITY["EPSG","54004"]]

2014-06-11 16:14:02,033 WARN [org.geoserver] - Ignoring store directory 'layergroups'

2014-06-11 16:14:02,595 WARN [org.geoserver] - Ignoring store directory 'styles'

2014-06-11 16:14:04,282 WARN [org.geoserver] - Error connecting to 'Taxlots'. Disabling.

2014-06-11 16:14:04,876 WARN [org.geoserver] - Ignoring store directory 'styles'

2014-06-11 16:14:07,439 WARN [org.geoserver] - Ignoring store directory 'layergroups'

2014-06-11 16:14:07,439 WARN [org.geoserver] - Ignoring store directory 'styles'

2014-06-11 16:14:07,501 ERROR [geotools.jdbc] - There's code using JDBC based datastore and not disposing them. This may lead to temporary loss of database connections. Please make sure all data access code calls DataStore.dispose() before freeing all references to it

2014-06-11 16:14:08,829 WARN [org.geoserver] - Ignoring store directory 'styles'

2014-06-11 16:14:09,454 WARN [org.geoserver] - Ignoring store directory 'styles'

2014-06-11 16:15:04,873 WARN [support.DisposableBeanAdapter] - Invocation of destroy
method failed on bean with name 'geoServerLoader': java.lang.NullPointerException

2014-06-11 16:15:04,873 ERROR [context.ContextLoader] - Context initialization failed
org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'gwcFacade' defined in URL [jar:file:/E:/path_to_install/OpenGeo%20Suite/webapps/geoserver/WEB-INF/lib/gwc-2.2-SNAPSHOT.jar!/applicationContext.xml]: Cannot resolve reference to bean 'gwcStorageBroker' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'gwcStorageBroker' defined in URL [jar:file:/E:/path_to_install/OpenGeo%20Suite/webapps/geoserver/WEB-INF/lib/gwc-2.2-SNAPSHOT.jar!/geowebcache-core-context.xml]: Cannot resolve reference to bean 'gwcMetaStore' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'gwcMetaStore' defined in URL [jar:file:/E:/path_to_install/OpenGeo%20Suite/webapps/geoserver/WEB-INF/lib/gwc-2.2-SNAPSHOT.jar!/geowebcache-core-context.xml]: Instantiation of bean failed; nested exception is org.springframework.beans.BeanInstantiationException: Could not instantiate bean class [org.geowebcache.storage.metastore.jdbc.JDBCMetaBackend]: Constructor threw exception; nested exception is java.lang.OutOfMemoryError: Java heap space


texts: How to add a layer, from ArcGIS Online, to Web AppBuilder at run-time

texts:
This is similar to my previous question How to add a layer to Web AppBuilder at run-time

But in this case I want to access layers from ArcGIS Online, similar to the way that the > Add > Search For Layers function works from the default ArcGIS Online map viewer:

This brings up a window listing layers stored in my ArcGIS Online account, and allows them to be added to the map.

Is a similar function available in Web AppBuilder?


texts: Does QGIS support Galaxy Tab 4 10.1 4G LTE?

texts: I need a tablet for outdoor work and i would like to know if this device "Samsung Galaxy Tab 4 10.1 4G LTE" works fine with QGIS mobile application.


texts: The general public and their ability to use GIS applications - do we still need paper maps?

texts:
I am endeavouring to convince others in my organisation that GIS tools, in particular digitising tools, CAN be used by the general public and we should promote their use instead of using paper maps.  This would save money (postal fees) and time (the time it takes to post maps) and I suppose it's better for environment (not using paper).

The counter argument is that people feel more at home with a pen/pencil and sketch out areas/polygons etc on a paper map, they also tend to be more accurate.  My business requires applicants to be fairly accurate with sketching out boundaries against a background map.

What do you think I can do to persuade those disbelievers? :) How can we improve sketching tools for non-GIS users so they make accurate sketches on maps.  What tips/techniques do you know of, have you faced this problem yourself??  Are we doomed to using paper maps for ever more!?  Your thoughts are most welcome.

PS - the GIS application is likely to be web based and open to the general public.  Probable tech would be ESRI Javascript API - with cached up, optomised map services, sensible scale levels and clear symbology etc.


texts: ArcGIS Online - Migrating items between accounts

texts: Any recommendation for moving the applications/items between ArcGIS Online accounts? I was trying to use ArcGIS Online Assistant, which uses the ArcGIS REST API to provide items migration, but it will break the links and content becomes lost by using this tool.


texts:
I'm very new to mint/linux/GRASS and have no clue how to deal with these errors.  Can I just change the URL somewhere?

Here is the simplest reproducible command to produce a failure:
    sudo apt-get install mysql-client-core-5.5

Here is the output from my attempt to run/install grass:

mint@mint ~ $ grass
The program 'grass' is currently not installed. You can install it by typing:
sudo apt-get install grass-core
mint@mint ~ $ sudo apt-get install grass-core
Reading package lists... Done
Building dependency tree
Reading state information... Done
The following extra packages will be installed:
  gdal-bin libarmadillo3 libdap11 libdapclient3 libepsilon0 libfreexl1
  libgdal1 libgeos-3.3.3 libgeos-c1 libgfortran3 libhdf4-0-alt libhdf5-7
  liblapack3 libmysqlclient18 libnetcdfc7 libodbc1 libogdi3.2 libpq5 libproj0
  libspatialite3 libxerces-c28 mysql-common odbcinst odbcinst1debian2 proj-bin
  proj-data python-numpy
Suggested packages:
  python-gdal grass-doc grass-gui grass-dev e00compr avce00 gnuplot gpsbabel
  gpstrans netpbm python-rpy2 xml2 libhdf4-doc libhdf4-alt-dev hdf4-tools
  libnetcdf4 libmyodbc odbc-postgresql tdsodbc unixodbc-bin ogdi-bin
  proj-ps-doc python-numpy-doc python-numpy-dbg python-nose python-dev
  gfortran
The following NEW packages will be installed:
  gdal-bin grass-core libarmadillo3 libdap11 libdapclient3 libepsilon0
  libfreexl1 libgdal1 libgeos-3.3.3 libgeos-c1 libgfortran3 libhdf4-0-alt
  libhdf5-7 liblapack3 libmysqlclient18 libnetcdfc7 libodbc1 libogdi3.2 libpq5
  libproj0 libspatialite3 libxerces-c28 mysql-common odbcinst odbcinst1debian2
  proj-bin proj-data python-numpy
0 upgraded, 28 newly installed, 0 to remove and 130 not upgraded.
Need to get 1,043 kB/33.6 MB of archives.
After this operation, 94.8 MB of additional disk space will be used.
Do you want to continue [Y/n]? Y
Err http://archive.ubuntu.com/ubuntu/ quantal-updates/main mysql-common all 5.5.28-0ubuntu0.12.10.1
  404  Not Found [IP: 91.189.92.156 80]
Err http://archive.ubuntu.com/ubuntu/ quantal-updates/main libpq5 amd64 9.1.6-1ubuntu1
  404  Not Found [IP: 91.189.92.156 80]
Err http://security.ubuntu.com/ubuntu/ quantal-security/main mysql-common all 5.5.28-0ubuntu0.12.10.1
  404  Not Found [IP: 91.189.92.200 80]
Err http://security.ubuntu.com/ubuntu/ quantal-security/main libmysqlclient18 amd64 5.5.28-0ubuntu0.12.10.1
  404  Not Found [IP: 91.189.92.200 80]
Failed to fetch http://security.ubuntu.com/ubuntu/pool/main/m/mysql-5.5/mysql-common_5.5.28-0ubuntu0.12.10.1_all.deb  404  Not Found [IP: 91.189.92.200 80]
Failed to fetch http://security.ubuntu.com/ubuntu/pool/main/m/mysql-5.5/libmysqlclient18_5.5.28-0ubuntu0.12.10.1_amd64.deb  404  Not Found [IP: 91.189.92.200 80]
Failed to fetch http://archive.ubuntu.com/ubuntu/pool/main/p/postgresql-9.1/libpq5_9.1.6-1ubuntu1_amd64.deb  404  Not Found [IP: 91.189.92.156 80]
E: Unable to fetch some archives, maybe run apt-get update or try with --fix-missing?
mint@mint ~ $


texts: what's on-the-fly processing in arcgis? saw this terminology while doing some research on mosaic dataset. Looking for a simple(less-complicated) answer.


texts: Intersecting vector data with raster data in PostGIS?

texts:
I have Open Street Map [OSM] data imported so that I can access the geometry of buildings.
I have also loaded SRTM digital height data from NASA from a GeoTIFF raster. I have projected the raster into EPSG:3857 (the same as OSM) from the original EPSG:4326.

For the SRTM tiles that I am interested in, I have warped them into 1 meter resolution with cubic spline interpolation via ST_Resample and produced smaller tiles from that with ST_Tile(as a single INSERT .. SELECT which is very handy).

Now I run into a problem:

The ST_Intersection() of the resampled tile and the building polygon simply produces a square bounding box, which is not what I want.  I want to get at the exact x,y points on the raster for the intersection because I'm using this data for generating 3d voxel scenes (1 pixel = one voxel) and I want to produce a 2d map of the rendered area to pixel-point precision as well.

What I am currently doing:

I use ST_Clip() to clip the raster by the geometry of the buildings, which puts NODATA data into the returned clipped raster everywhere except where the building is actually situated.  This is like using an template to "burn" the image of the building into the raster, leaving the points inside the building completely intact while everything else is burned away.

I dump the resulting clipped raster with ST_PixelAsPoints() or ST_PixelAsCentroids() which gives me x,y positions on the raster where the building is located.

Is this the right approach?

I have heard of something called a 'pointcloud' but PostGIS doesn't seem to have this available.  Converting the raster into polygons doesn't work, because the x,y data from the raster is lost.


texts: How to choose scale endpoints for` tm_fill( )`?

texts:
I'm trying to map a continuous scale with tm_fill in tmap:

tm_shape(shape)+tm_fill(col="value",pallete="seq",contrast=c(0,1)+tm_facets(by=group)

I want to make several such maps for different years but I want to have the scale for the fill variable to be between 0 and 1 for all of them. I thought that contrast would do that, but every map has a scale that ranges between minimum and maximum value.


texts:
I have a pretty simple script which gets 2 parameters from the user:

1)Geodatabase location

2)Spatial reference

It then creates a dataset with 6 feature classes. It also adds a large amount of fields to each each feature class, and by large I mean no more than 20.

It also adds defaults to certain fields.

Could I speed up the running of the tool by organizing the code differently?

Right now I have it so it creates the dataset, then creates the first feature class and then the fields it needs and any defaults for those fields.

Could it be quicker to have it make all the features classes at once, then all the fields, then all the defaults? Or would it make no difference?

Currently, it's taking about 4 minutes to do everything.


texts: Get the probable right location from multiple points

texts:
I have a location collection tool which send me location data of our survey. Surveyors send multiple data of a single item. Such as they send numarous data for a school like the following image. From the image we can clearly identify few locations those are actually not collected properly or inaccurate location came due to GPS accuracy. Here I have marked red circle which location need to be excluded and the green circle indicates this is the probable location. IF we create a centroid based on all points then it will not indicates the probable correct location due to few highly inaccurate data. I wonder, if any method or tool can exclude the wrong location and create a cluster based on densed location points to create a centroid then the centroid can be established as probable right location of the surveyed item. I have millions of item those have hundreds of location each. I came to know that k-means do something like that but it does not do on a single item. I have four column in table.

fid
item_id
item_name
geom

Is there any method to create item wise clustering where every cluster have a value which indicates its containing points to get the probable right location?


texts: Append not merge multipolygon to existing multipolygon in QGIS from PostGIS database

texts:
We have areas (merged polygons) with an unique id in a PostGIS database. The unique id (serial) is referenced by other datasets (actions that happen in these areas). Sometimes we need to append a few new polygons to the existing areas. When we try to merge them simply together, this doesn't work because the original area (with id) will be deleted to create a new one, but this isn't possible when the old area has actions that rely on the id.
Is there no other possibility to append the new polygons to the old one without changing (deleting) the original id? Just like updating the geometry of the dataset.


texts: What algorithm does QGIS use to calculate Nearest Neighbor Index?

texts:
I'm carrying out nearest neighbor analysis on a set of geological features in QGIS, trying to find out whether they're clustered or not. I've calculated the nearest neighbor distances using the 'Distance Matrix' tool and I'm now exploring the 'Nearest Neighbor Analysis' tool.

Does anyone out there know what the algorithms behind the NN Analysis tool are? What method does it use to generate its expected distance and NN index (comparison with Poisson points?)? How does it measure the area? e.g. does it use minimum convex hull or minimum bounding rectangle?

I became a bit suspicious when I tried the tool on a set of randomly generated points, it returned NN index values of 0.81 and a Z-score of -2.5, indicating clustering in what should be randomly distributed points.

Can anyone shed any light on this? I've tried to find more information in the QGIS documentation, but there's not information about the algorithms that are used to calculate these values.


texts:
I have a large dataset which I need to display in an online map. We use Geocortex and ArcGIS Online services. I have heard of tile cache, but I can't work out what I need to do to create it.

Can I do it in ArcGIS Pro?


texts:
When I run a batch script using a GDAL/OGR program, I can detect if an error occured by checking ERRORLEVEL for the value 1.
When using GDAL/OGR in python or .NET, I get an exception when an error occurs.

Question:
How do I detect if a Warning was issued during the execution?

Why?
Even though a process can execute to the end despite of warnings, the output might be corrupt with respect to the data I expect. You could argue that I ought to check for these corruptions, but since I have control of the entire process, I would like to be able to stop when a warning is issued.


texts: QGIS 3 Plugin Cancel Edits is crashing plugin

texts:
I was using, this method in QGIS2.

self.iface.actionCancelEdits().trigger()

This method was used, on canceling the edits in the iface. But this method is crashing my plugin on QGIS 3. What is the alternative to this method or am i doing somethign wrong with this method ? I am also using below syntax in QGIS 3, but it is also not working.

self.dlg.showNormal()


texts: PyQGIS update field in line dataset when it is within a polygon

texts:
I am trying to update a field in a (Road centerline) table based on a value from another table when it is within a Polygon (Locality boundaries).

Based on posts like QGIS layer.geometry.intersection() not finding intersections between layers I have the following code. Both datasets are in EPSG:28355. It never prints ntext=b['Locality_N'] - so obviously something is wrong with the if within statement.

    import re

    def Update_Field(fieldName, ntext): # Function to simplify updating fields
        print (ntext)
        fieldIndex = layer.fields().indexFromName(fieldName)
        dpr.changeAttributeValues({f.id(): {fieldIndex: ntext}})
        layer.commitChanges()

    RoadType_dict = {'Rd':'Road', 'Avenue': 'Ave'}
    layer = iface.activeLayer()
    selection = layer.getFeatures() #all features
    #selection = layer.selectedFeatures() #only selected features

    #Locality Data
    Locality_layer=QgsVectorLayer(r'Z:\My Drive\Mangoesmapping\Spatial Projects\Gen_Library\DSC\vector\DSC_Localities_areas.shp', "Locality", "ogr")

    dpr = layer.dataProvider()

    for f in selection:
        print (f.id())
        text=f['RoadMntnc']
        print (text)
        stext=re.split('_+', text)

    ##Update Road Type
    #    ntext=RoadType_dict.get(stext[-1:][0])
    #    fieldName='Type'
    #    Update_Field(fieldName,ntext)

    ##Update LocalityName
        for b in Locality_layer.getFeatures():
            if f.geometry().within(b.geometry()):
                print (b.geometry())
                ntext=b['Locality_N']
                fieldName="Locality"
                Update_Field(fieldName,ntext)

    layer.commitChanges()
    layer.updateFields()

I have commented out the code that works but is not required for these tests.

--- Based on another answer I also tried

e = QgsExpression( "geomwithin('Locality_layer','Locality_N') ")
f[fieldName] = e.evaluate()
layer.updateFeature(f)

But the values aren't updated in the table.


texts: Improve database performance for geospatial

texts:
I have a database with osm data imported with osmosis tool. I want to search for the nearest POIs from a given point. For this purpose, I perform this query:

select id from nodes
where geom='0101000020E6100000546C713216920DC065767176C62C4540'

After this I would like to search tags for the selected node in node_tags table but my problem is that the database has 31 million rows, and the query takes more than 30 seconds. I have this columns in node tables:

id, version, user_id, tstamp, changeset_id, geom

I run the following sentence to ensure index in "geom" column:

CREATE INDEX idx_nodes_geom
ON nodes USING gist (geom );

But it doesn't work. Query takes more than 30 seconds.

I have the following configuration in PostgreSQL database:

Shared buffers: 768 MB
Work_mem: 32MB
maintenance_work_mem: 256 MB
wal_buffers: 2 MB
checkpoint_segments: 12
random_page_cost: 2.0
seq_page_cost: 1.0

Any idea to improve it?


texts: Using OpenLayers (with OpenStreetMaps) and PostGIS for rendering data acquired from GPS tracking device

texts:
I am researching and planning to acquire data from GPS device and store them in PostGIS database, then use geodjango to process/work them out and finally chart them on a map with OpenLayers.
Will these technologies be enough for my plan, since I don't want my maps?
I want to use OpenStreetMap and plot my points on an OpenLayers' vector layer.
I also want to implement real time tracking on a map, I'll query the PostGIS database in real-time using AJAX and then recreate objects on a map. I just cannot find any good examples to do this with OpenLayers. Can I use Jquery for this and combine it with OpenLayers since I'm good at using Jquery's AJAX methods.
I am still researching these technologies, but OpenLayers/geodjango/postGIS seems like a great combo for my needs. I found literature for almost everything that I'm interested in these 3, only that I lack the examples of real time charting and refreshing objects on OpenLayers' vector layer.


texts: Geoserver raster as wms has blue-ish appearance when black background is visible

texts:
I have a problem with a blue-ish overlay color over my WMS layer that only is visible if there is part of the image background (nodata area) in the requested map extent.

The original image is a facsimile of an old map. It is georeferenced, resulting in approx. 45° rotation. Thus, the GeoTiff has large "no data" areas in the corners.

See original image, which is displayed correctly by geoserver's wms if there is no nodata visible in the requested extent:

...and if the nodata area on the bottom left is inside the viewport, it looks like this:

The SLD file used is

<?xml version="1.0" encoding="UTF-8"?><sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
  <sld:NamedLayer>
    <sld:Name>B0448_7c</sld:Name>
    <sld:UserStyle>
      <sld:Name>B0448_7c</sld:Name>
      <sld:FeatureTypeStyle>
        <sld:Rule>
          <sld:RasterSymbolizer>
            <sld:Geometry>
              <ogc:PropertyName>grid</ogc:PropertyName>
            </sld:Geometry>
            <sld:ChannelSelection>
              <sld:RedChannel>
                <sld:SourceChannelName>1</sld:SourceChannelName>
                <sld:ContrastEnhancement>
                  <sld:GammaValue>1.0</sld:GammaValue>
                </sld:ContrastEnhancement>
              </sld:RedChannel>
              <sld:GreenChannel>
                <sld:SourceChannelName>2</sld:SourceChannelName>
                <sld:ContrastEnhancement>
                  <sld:GammaValue>1.0</sld:GammaValue>
                </sld:ContrastEnhancement>
              </sld:GreenChannel>
              <sld:BlueChannel>
                <sld:SourceChannelName>3</sld:SourceChannelName>
                <sld:ContrastEnhancement>
                  <sld:GammaValue>1.0</sld:GammaValue>
                </sld:ContrastEnhancement>
              </sld:BlueChannel>
            </sld:ChannelSelection>
            <sld:ContrastEnhancement/>
          </sld:RasterSymbolizer>
        </sld:Rule>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:NamedLayer>
</sld:StyledLayerDescriptor>

If I set the layers InputTransparentColor to #000000 the wms will not show anything (only a white background in OpenLayers preview).

Unfortunately I cannot share the geotiff because it contains proprietary data.

How can I get rid of the blue?

UPDATE: I have now narrowed my problem down to (I think) a custom color stretch applied for every rendered tile and not for the whole raster. When using the options toolbar in the openlayers preview and viewing as tiled, all tiles will have a slightly different color scheme. Inner tiles (which dont touch the nodata areas on the edges) are near the original color, tiles containing the edges are blue-ish.


texts: Labeling layers with Python in QGIS. LinePlacementFlags

texts:
How could I set LinePlacementFlags with Python in QGIS?
My code is something like:

palyr = QgsPalLayerSettings()
palyr.readFromLayer(layer)
palyr.enabled = True
palyr.fieldName = "attribute"

Then I've tried:

palyr.placement= QgsPalLayerSettings.Line

And my labels are on the line. I want them above the line, but I can't figure how to do this.
I found there is LinePlacementFlags type. But

palyr.LinePlacementFlags= QgsPalLayerSettings.AboveLine

didn't work.


texts: RuntimeError: cannot open  from ArcPy search cursor?

texts:
import arcpy
from arcpy import env
env.workspace= "C:/Esripress/Python/Data/Excercise08"
fc = "rivers.shp"
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@LENGTH"])
length =0
for row in cursor:
    length = length + row[0]
print length

I am getting the following error:

"Traceback (most recent call last):
File "C:\EsriPress\Python\Data\Exercise08\Results\geometry.py", line 5, in <module>
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@LENGTH"])
RuntimeError: cannot open 'rivers.shp'"

My data is on C:\EsriPress\Python\Data\Exercise08\Results


texts:
I am using gdal_merge to create an raster from around 250 files (rasters). when I run the command i got gdal_merge was unexpected at this time

The command is correct as if I split the files it works fine

gdal_merge -o "Z:\30cm Imagery 2015 original\test.tif" "Z:\50cm Imagery 2015 original\jp2\15003800.jp2" "Z:\50cm Imagery 2015 original\jp2\15003850.jp2"


texts: How to use gdal.openEx with flags and parameters?

texts:
I'm new in python. I have a class which opens SXF format, and i want the last string

 self.sxf_ds = gdal.OpenEx(self.sxf_path)

to be opened with flag 'GDAL_OF_VERBOSE_ERROR' but the  nOpenFlags is integer not the string. Where i can find the NUMBER of this flag??.
 In gdal.py for openEx function there is a comment:
"OpenEx(char const * utf8_path, unsigned int nOpenFlags=0, char ** allowed_drivers=None, char ** open_options=None, char ** sibling_files=None) -> Dataset"

In gdal help only definition is written:
nOpenFlags a combination of GDAL_OF_ flags that may be combined through logical or operator.

Driver kind:GDAL_OF_RASTER for raster drivers, GDAL_OF_VECTOR for vector drivers. If none of the value is specified, both kinds are implied.GDAL_OF_RASTER for raster drivers, GDAL_OF_VECTOR for vector drivers. If none of the value is specified, both kinds are implied.
Access mode: GDAL_OF_READONLY (exclusive)or GDAL_OF_UPDATE.
Shared mode: GDAL_OF_SHARED. If set, it allows the sharing of GDALDataset handles for a dataset with other callers that have set GDAL_OF_SHARED. In particular, GDALOpenEx() will first consult its list of currently open and shared GDALDataset's, and if the GetDescription() name for one exactly matches the pszFilename passed to GDALOpenEx() it will be referenced and returned, if GDALOpenEx() is called from the same thread.
Verbose error: GDAL_OF_VERBOSE_ERROR. If set, a failed attempt to open the file will lead to an error message to be reported.

here is a link
to gdal docs.

If i write self.sxf_ds = gdal.OpenEx(self.sxf_path,'GDAL_OF_VERBOSE_ERROR')
i'm getting an obvious error: TypeError: in method 'OpenEx', argument 2 of type 'unsigned int'. I suppose i need something like that self.sxf_ds = gdal.OpenEx(self.sxf_path,1)or self.sxf_ds = gdal.OpenEx(self.sxf_path,2)

class sxf2dxf(object):
    """docstring for sxf2dxf"""

    def __init__(self, sxf_path, rsc_path, dxf_style_path, dxf_out_path):
        self.sxf_path = sxf_path
        self.rsc_path = rsc_path
        self.dxf_style_path = dxf_style_path
        self.dxf_out_path = dxf_out_path

        self.cls = RSC.parse(rsc_path, cls_type='geom_type')
        tools.MyPrint().pprint(self.cls)
        gdal.SetConfigOption('SXF_RSC_FILENAME', self.rsc_path.encode('utf-8'))
        # gdal.SetConfigOption('SXF_LAYER_FULLNAME', False)  # usage of long layer names

        self.sxf_ds = gdal.OpenEx(self.sxf_path)


texts:
I am trying to add a trail line to a story map. I can not find where in the story map it allows for adding a layer.

Anyone have suggestions?


texts:
My current pet project involves creating a web service that uses the CGIAR SRTM worldwide elevation data set. Each request to the service will query for elevation data from a small geographic area.

Storing all the data in fashion that can be queried fast turns out to be difficult. My best idea is to keep the data in GeoTIFF files and create an index linking lat/lon coordinates to files. Which leads to my current main problem:

How do I get the elevation values out of a GeoTIFF file? I know the projection, scale and orientation of the GeoTIFF files (i.e. I know what lat/lon values each pixel represents). How can I figure out what elevation that pixel refers to? I have to use PHP and thought I could use the ImageMagick function for reading a pixel's color. Unfortunately, it returns R=2 G=2 B=2 for all pixels.

update

To answer someone else's question below, here is my solution.

I went with the library by Bob Osola suggested in the accepted answer below and made the following changes to adjust it to SRTMv4.1 file format (the one available from the CGIAR website). Update the constants to:

const NUM_DATA_ROWS = 6000;     // the number of data rows in the file ( = ImageLength tag value)
const NUM_DATA_COLS = 6000;     // the number of data columns in the file ( = ImageWidth tag value)
const DEGREES_PER_TILE = 5;     // each tile is 5 x 5 degrees of lat/lon
const PIXEL_DIST = 0.000833333; // the distance represented by one pixel (0 degrees 0 mins 3 secs of arc = 1/1200)
const STRIPOFFSETS = 0x44bbd5c; // the offset address of the 'StripOffsets' tag

In addition there was a bug in the computation of a longitude. To fix it, replace

$topleftLon = 180 - (($tileRefHoriz -1) * self::DEGREES_PER_TILE);
if ($lon < 0) {
    $topleftLon = -$topleftLon;
}

by

$topleftLon = (($tileRefHoriz -1) * self::DEGREES_PER_TILE) - 180;

That should be it.


texts: GeoServer WFS 2.0 GetFeature query with wildcard for typeNames/featureTypes

texts:
Just looking through GeoServer docs trying to find a way to construct a WFS 2.0 query to get all features within a BBOX - regardless of feature type.

The service I'm consuming will be changing over time, so I'd prefer to use a wildcard type argument to cover all feature types in the request rather than hard-code them..

Hopefully without having to query the service prior just to determine the feature type params.

Have tried the usual syntax, e.g. typeName=* but no good.

Is there something I'm missing?


texts:
I create a RelQueryTable like this to join a GeoFeatureLayer to an external table.  Everything seems to work (no errors, I can look in the table and see the correct fields.) Silly question: what do I return from this that I can use as a GeoFeatureLayer? When I try to cast the table to GeoFeatureLayer it will not work. It was hard enough to riddle out how this even worked (and I still have many unanswered questions that the bathtub documentation doesn't help with while I'm adrift at sea.)

Does anyone have intelligent samples for these APIs? (Instead of changing 9.x to 10.x they should have released some code.)

Type memoryRelationshipClassFactory = Type.GetTypeFromProgID
    ( "esriGeodatabase.MemoryRelationshipClassFactory" );
IMemoryRelationshipClassFactory memoryRelationshipFactory =
    Activator.CreateInstance( memoryRelationshipClassFactory ) as
    IMemoryRelationshipClassFactory;

IRelationshipClass relationshipClass = memoryRelationshipFactory.Open(
     "Join",
geoFeatureLayer.DisplayFeatureClass as IObjectClass,
layerField,
table as IObjectClass,
tableField,
"forward",
"backward",
esriRelCardinality.esriRelCardinalityOneToOne );

Type rqtFactoryType = Type.GetTypeFromProgID( "esriGeodatabase.RelQueryTableFactory" );
IRelQueryTableFactory rqtFactory = Activator.CreateInstance( rqtFactoryType ) as
    IRelQueryTableFactory;

IRelQueryTable relQueryTable = rqtFactory.Open(
    relationshipClass,
    true,
    null,
    null,
    EXTERNAL_FIELDS,
    true,
    false ) as IRelQueryTable;


texts: Merge CSV file to Vector Layer using common attributes in QGIS python scripting

texts:
I have a vector layer and a delimited text (csv) file -- they both have the common attribute attribid -- how can I add the columns of the csv file to the vector layer merging on the the attribid key using a python script?

I'm aware that I can use MMQGIS's combine feature to do this, but I could not get it to work using a python script.


texts:
I've imported a shapefile into R, and joined it to a table. My shapefile contains  all the census ids, while my table only contains selected census ids. I'm now trying delete all the rows didn't get a match.

This is what my dataset looks like (I'm trying to remove all the rows with NA, so the last two would need to be removed):

  CTUID   Cluster   Average
5350007.01  1       124.53
5350007.02  1       234.87
5350010.01  4       110.11
5350010.02  5       187.68
5350001     NA        NA
5350002     NA        NA

I've tried using this line code:

shape2[!(rowSums(is.na(shape2))==NCOL(shape2)),]

Which gave me this error:

Error in rowSums(is.na(shape2)) :
  'x' must be an array of at least two dimensions
In addition: Warning message:
In is.na(shape2) : is.na() applied to non-(list or vector) of type 'S4'

I'm not very proficient in R, so any help would be really appreciated. If you could include a brief explanation that would be fantastic.


texts: Problem running intersect vector processing with line and polygon shapefiles

texts: I have been trying to run intersect on two shapefiles, one is line, the other is polygon. Although processing appears to take place and a new shapefile is created and added, nothing is visible. I have checked their projections which are set the same.


texts:
I have to draw in ArcGis a graph which has the end points of each segment with a circle. The circle has to be unique at the same point when two or more end points are coincident.

I’m using the font of ESRI Default Marker (Unicode 80) and the circle is not centred at the end point, but it is shifted of the ray's length along the segment. So, the result is that I have two different circles for each segment.
If I give an offset at the position of the symbol the result is still two or more slightly different circles, but not coincident.

Is there a way to obtain an unique circle when more segments are coincident at the same end point?

Thanks for any help! I’m using ArcGis 9.3


texts: How to access a feature by specific field value using python ogr in an efficient way?

texts:
I have one point layer and one polygon layer with a common identifier attribute field called "N", same objects amounts but Feature IDs (FID) are different.

I have to iterate on the first one and what I want is to get the FID from the second using an attribute join (access to FID from the polygon using the specific field value of "N" as a bridge between both layers.

for feature in pointsLayer
    bridge = feature.GetField('N')

    # Here I want to select my feature from my other layer
    # which has the same attribute value "N"
    polygonsLayer.GetFeature(...)

As I read for now, I have 2 solutions, SQL request or create a matching table N-idPolygon?

Any hint what would be more efficient? We talk about 10000 features max in each layer.


texts: Im using leaflet, Is there a plugin that allows it to print more than one map in the same pdf like this example .


texts: Clipping multiple resampled rasters to same size using ArcGIS Desktop?

texts:
I have various rasters of environmental variables that were imported as GeoTiffs or NetCDFs. The rasters had varying cell size which I fixed by resampling them using the Snap Raster function in Environments > Output Extent to ensure that the cell borders of all rasters match up.

So far so good.

However, when I now try to clip all rasters to a rectangluar polygon some of the clipped rasters have an additional row or column, i.e. the clipped rasters have different row or colum sizes.

How do I clip all rasters to the same size using ArcGIS Desktop 10.3?


texts:
I'm using Environment Agency DTM LiDAR data which I import into Mathematica in ASCII format to produce ground surface models, sections, etc.

The process works fine, but I'm unclear whether the cell elevation values always relate to the easting and northing at the cell centroids or whether their plan positions are dictated by the header data XLLCORNER/XLLCENTER and YLLCORNER/YLLCENTER?


texts: Creating domains from specific values in table using ArcPy?

texts:
I have an .sde table that I'm trying to extract specific values from and create domains with them. The problem I'm running into is that there are multiple duplicate values in the table (which is fine, it's supposed to be like that due to different categories) and I'm only interested in the unique ones. I've created a dictionary valueDict which contains my domain names as the value and the application they're located in as the key. I then need to iterate over the values I'd like as my domains and insert them into their respective tables (the value from above).

Each application has a specific set of values that I'd like as domains. I'm also filtering the table to receive yet another specific view. I'm only printing the values below to see what I'm returning just for testing purposes.

In the example below, I'm getting back every value, including the duplicates.

...
InTable = "MyTable"
Fields = "DESCRIPTION", "TIME"

for key, value in valueDict.items():
    exp = "APPNAME = " + key + " and APR = 'EVENT'"
    with arcpy.da.SearchCursor(InTable, Fields, exp) as cursor:
        for i in cursor:
            print(key, value, i[0], i[1])


texts: Can joined tables be served by QGIS Server WFS?

texts:
I'm trying to figure out whether it's just me, or if I'm forgetting a setting. I have a shapefile which I make available for WFS-requests. Now I added a join (to a MySQL table).

I tried playing around a bit with permissions, and added all (insert/update/delete/select) for both the layer and the table. I also checked  "fields" in the layer - I see the fields I join there. All the WFS-fields are checked, so I am expecting all the fields to be available/listed in the WFS response.

But the joined ones don't show in in the WFS response.

QGis 2.6 on Debian.


texts: Using Rasterio, Matplotlib, and Numpy trying to calculate and plot NDVI of a raster image getting error: IndexError: list index out of range?

texts:
First I create two variables and store my raster images in them and I use MatPlotLib to plot a subset of the original "nirband" raster image. That works fine. I am getting an IndexError: list index out of range error now. I am basing my code off the following tutorial found on GitHub:

https://github.com/geohackweek/tutorial_contents/blob/ba5e9443137a9aca87cdcdcd70e9e6a237cc64ba/raster/notebooks/rasterio-landsat-aws.ipynb

Here is the first section of the code which plots fine without error:

import rasterio
import matplotlib.pyplot as plt
import numpy as np

nirband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF"

redband =r"LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF"

#rasterio.windows.Window(col_off, row_off, width, height)
window = rasterio.windows.Window(1024, 1024, 800, 600)

with rasterio.open(nirband) as src:
    subset = src.read(1, window=window)

plt.figure(figsize=(15,8))
plt.imshow(subset)
plt.title(f'Band 5 Subset\n{window}')
plt.xlabel('Column #')
plt.ylabel('Row #')

All the code before this point runs fine, I was able to plot the above code using MatPlotLib just fine. It's the following code which produces an error message where I am trying to compute and display the NDVI of that raster subset image:

with rasterio.open(redband) as src:
    profile = src.profile
    oviews = src.overviews(1) # list of overviews from biggest to smallest
    oview = oviews[1]  # Use second-highest resolution overview
    print('Decimation factor= {}'.format(oview))
    red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))

plt.imshow(red)
plt.colorbar()
plt.title('{}\nRed {}'.format(redband, red.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')

with rasterio.open(nirband) as src:
    oviews = src.overviews(1) # list of overviews from biggest to smallest
    oview = oviews[1]  # Use second-highest resolution overview
    nir = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))

plt.imshow(nir)
plt.colorbar()
plt.title('{}\nNIR {}'.format(nirband, nir.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')

def calc_ndvi(nir,red):
    '''Calculate NDVI from integer arrays'''
    nir = nir.astype('f4')
    red = red.astype('f4')
    ndvi = (nir - red) / (nir + red)
    return ndvi

ndvi = calc_ndvi(nir,red)
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar()
# maybe edit it
plt.title('NDVI')
plt.xlabel('Column #')
plt.ylabel('Row #')

Here is the complete traceback with error message:

runfile('C:/Users/new/Desktop/RasterNDVI.py', wdir='C:/Users/new/Desktop')
Traceback (most recent call last):

  File "<ipython-input-15-5d04fa0ce75f>", line 1, in <module>
    runfile('C:/Users/new/Desktop/RasterNDVI.py', wdir='C:/Users/new/Desktop')

  File "C:\Users\new\Miniconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile
    execfile(filename, namespace)

  File "C:\Users\new\Miniconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/new/Desktop/RasterNDVI.py", line 27, in <module>
    oview = oviews[1]  # Use second-highest resolution overview

IndexError: list index out of range

Attempt at saving NDVI as a GeoTIFF file:

with rasterio.open("LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF") as src:
    naip_data_ras = src.read()
    naip_meta = src.profile

with rasterio.open('MyExample.tif', 'w',**naip_meta) as dst:
    dst.write(naip_ndvi, window=window)


texts: Make multipoly only if two field values overlap?

texts:
So I want to combine all the polygons that share two different attribute values and don't really know where to start.

The selection shows you my problem, they both have the same classification (890) and belong to the same block (third column). I want these (and all other similar occurences) to be considered a single multipart polygon. How?

The regular singlepart to multipart tool won't do me any good since I will loose either classification field or the block id field.


texts: Error Trying to Run R Summary Statistics in QGIS 2.12.1

texts:
When I try to run the R script "Summary statistics" in QGIS 2.12.1 on a Windows 7/64 bit machine, I get the error message "error in if (!file.info(DSN)$isdir) DSN missing value where TRUE/FALSE needed," and the output file is empty.

Having read a comment that R doesn't work properly if installed in the c:\Program Files\ directory, I tried installing it in its own directory in c:, but that didn't make any difference. I've tried installing R with and without the 32 bit files, and with just the 32 bit files (no 64 bit files), but the result is always the same: no results. The same happens with the R Histogram script.

QGIS's native Basic Statistics routine does produce results. Can anyone help me get R scripts working in QGIS?


texts: Modifying raster values in a zone with a math expression

texts:
I'm currently doing an SIG project with the goal of calculating groundwater volumes in chalk and silex residuals aquifers.

I have point data of underground water, chalk roof, and floor elevations, as well as RS thickness. The problem is that the chalk floor data points are visibly taken from an interpolation while roof and water are taken from measurements. As a result there're some locations were the floor elevation exceeds roof and/or water ones.

I have made TINs from the data points and then created rasters from those. I would like to update the raster cells in the zone where floor elevation are not correct. Those zones are high slopes (either valley or tectonic features), so basically the problem is that floor elevations do not follow the slopes that roof and water do follow.

How could I update those floor cell values according a the profile/slope of water/roof values? I thought about applying a math expression but couldn't find any way to do so.


texts:
I'm using gdaldem color-relief to read a Uint32 raster and write an RGB tif. The input color_text_file has one line:

102216130    23:105:102

My input file is a single (large) pixel of value 102216130, matching the one line in the color text file.

If I run gdaldem with the -exact_color_entry flag, the output is a pixel with RGB 0 0 0. If I run it without the flag, the output is correct, a pixel with 23 105 102.

Why is this? It certainly doesn't matter for this raster (-exact_color_entry appears to be unnecessary), but for rasters with more lookup values in the color text file (700, in one case), the general case (i.e. not using -exact_color_entry) is important, as output colors do not exactly match my input.

The help page for gdaldem states that -exact_color_entry:

  use strict matching when searching in the color configuration file. If none matching color entry is found, the "0,0,0,0" RGBA quadruplet will be used

Can I edit this color configuration file? The colors that are being approximated seem normal (R 21, G 191, B 217 for example) but the exact values are not being written to tiff.

Alternatively, is there a better tool than gdaldem for essentially reclassifying one raster to three bands?


texts: Converting .obj 3D files to .tiff raster (or reading .obj on QGIS)

texts:
I received a 3d terrain model captured with a drone in .obj format. I am trying to open it in QGIS to use the elevation information as a georeferenced raster, but so far I can't open it or transform it to GeoTIFF. I can visualize the model superimposed with an orthophoto using Acute 3D viewer, but this program doesn't have the option to export to raster.
I have found multiple answers online about how to convert a DEM or TIFF to .obj, but not the other way around. I prefer to ask here before downloading and installing multiple software to find one that does the job.
There is a similar question unanswered, but it doesn't include the ability to transform .obj to raster using other software different from QGIS.
So, bottom line, is there a way to transform .obj to GeoTIFF? or else, is there a way to open an .obj file in QGIS?


texts: Unable to view WFS as Vector Layer (OpenLayers 3)

texts:
I am trying to add a WFS to map using OpenLayers 3. There is no error message and no layer displayed on map. GeoJSON object is received in response.  Styling object is also getting initialized. Need help to understand if I am missing anything.

<!doctype html>

<link rel="stylesheet"
 href="http://openlayers.org/en/v3.10.1/css/ol.css" type="text/css">

<title>MAP</title>

Sample Map


texts: Does Leaflet always work in decimal degrees?

texts:
I'm looking into Leaflet and I'm surprised that it doesn't (seem to) specify coordinates in Web Mercator like other popular APIs, but instead uses decimal degrees.

Are all coordinates in Leaflet specified in decimal degrees? Is there a function to translate into other coordinate systems (as in the geographicToWebMercator option in ArcGIS)?

Thanks


texts: convert .shp to .gpx with selected attributes

texts:
I am working with QGIS 2.6/2.8 (1.8 is installed too).

Right now I am planning an inventory and want to convert a .shp file to a .gpx including the Column "fid" from the attribute table. Is there a way to do it? Or can I copy it later on?

There a more columns in the table of my .shp but I only need the "fid".


texts:
Im unable to get EucAllocation function working in ArcMap 10.5.1 or via Arcpy. I've used user guides for info but with no joy - http://desktop.arcgis.com/en/arcmap/10.5/tools/spatial-analyst-toolbox/euclidean-allocation.htm. I have noticed that using the "copy code snippet" from ArcMap is different to the format in above link.

code snippet (using defaults) = arcpy.gp.EucAllocation_sa(input,output, "", "", "5", "Value", "", "")

I have added the function where problematic part of the python script is. The code is generating a floodplain dataset from a DEM raster.

def create_floodplain(root_folder, in_gdb_name, working_gdb_name):

    # use output of get_raster function as input
    raster = get_raster(root_folder, in_gdb_name)

    working_gdb = os.path.join(root_folder, working_gdb_name)

    output_root = os.path.join(working_gdb, raster)
    #1 Raster Calc - convert DEM from floating point to int
    output_p1 = os.path.join(output_root + "_p1_rastercalc_TEST")

    print "Step 1. Running Raster Calculator...."
    rast_calc = Int(str(raster))
    # save output
    rast_calc.save(output_p1)
    print "Raster Calculator complete - %s" % (os.path.basename(output_p1))

    #2 Extract by Mask - create water boundary (ridgeline) raster
    input_p2 = output_p1
    mask_fc = os.path.join(working_gdb, "WATER_POLY_Clip16")
    output_p2 = os.path.join(output_root + "_p2_extractmask_TEST")

    print "Step 2. Running Extract By Mask...."
    ext_by_mask = ExtractByMask(input_p2, mask_fc)
    # save output
    ext_by_mask.save(output_p2)
    print "Extract By Mask complete - %s" % (os.path.basename(output_p2))

    #3 Euclidean Allocation - calculates, for each cell, the elevation value of the nearest ridgeline
    input_p3 = output_p2
    mask_fc = os.path.join(working_gdb, "WATER_POLY_Clip16")
    output_p3 = os.path.join(output_root + "_p3_euclideanallocation_TEST")

    print "Step 3. Running Euclidean Allocation...."
    euc_allocate = EucAllocation(input_p3, "", "", "5", "Value", "", "")
    # save output
    eucAllocate.save(output_p3)
    print "Euclidean Allocation complete - %s" % (os.path.basename(output_p3))

    #4 Raster Calc - create vertical buffer where 3 is the required buffer height (m)
    input_p4 = output_p3
    output_p4 = os.path.join(output_root + "_p4_rastercalc_TEST")

    print "Step 4. Running Raster Calculator...."
    rast_calc = raster >= (input_p4 + 3)
    # save output
    rast_calc.save(output_p4)
    print "Raster Calculator complete - %s" % (os.path.basename(output_p4))

The running of the EucAllocation (Step 3) in the above code is causing the issue as it is resulting in the below error:

  ExecuteError: ERROR 999999: Error executing function. A column was
  specified that does not exist. A column was specified that does not
  exist. The table was not found. [VAT_EucAllo_Suir1] No spatial
  reference exists. ERROR 010296: Error in writing raster
  C:\Users\XXX\AppData\Local\Temp\t_t83.tif. Distance mapping is failed.
  ERROR 010067: Error in executing grid expression.
  Failed to execute(EucAllocation).


texts:
I am new to ArcGIS and geoids/ellipsoids, and I want to convert a dataset in EGM96 (the USGS SRTM DEMs) to NAVD88. I have the SRTM dataset from NASA downloaded and open in ArcGIS as .hgt files.

I have downloaded VDatum from NOAA, which lists the EGM96 geoid on its webpage, but when I try to use the tool there are no EGM geoids available to select.

I have also found this https://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/binary/binarygeoid.html  which seems to be code written in FORTRAN, available to the public, to allow conversion of EGM96 to NAVD88. However the input doesn't seem to be .hgt files, which is what the SRTM dataset is. Also, this looks difficult to implement (and I'm not too familiar with FORTRAN77) and I'm wondering if there is something more recent.

Does anyone have a (relatively) simple approach to this? I don't think ArcGIS will perform the conversion. I have version 10.6. I haven't tried it myself as I'm new to ArcGIS but this is what I've read so far.

I am new to this topic so writing my own code for this is not something I can try at the moment (would take months).


texts: I need to extrude the 3D Borewell data using the inclined angle in ArcScene. I've seen only the options for extruding them in one direction(i.e., perpendicular downwards  or upwards). I got point data with Direction, distance and I need to use both the attributes for extruding below the surface.


texts: PostGIS: Extracting borders between contiguous (but not perfectly aligned) polygons

texts:
I have a (time-varying) shapefile with country polygons from which I want to extract borders between neighboring countries as a separate shapefile, using PostGIS. I managed to do this for most countries, but it did not work for all of them:

I did the following:

Transformed the polygons to lines:

CREATE TABLE t.test AS SELECT a.id, a.name, a."start", a."end", ST_boundary(wkb_geometry) AS geom from t.testfile;

Extracted the shared borders:

CREATE TABLE t.test2 AS SELECT a.id, a."name",a."start", a."end", b.id, b."name", b."start", b."end", ST_Intersection(a.geom, b.geom)
FROM t.test AS a, t.test AS b WHERE ST_Intersects(a.geom, b.geom) AND a.id != b.id AND (a.start, a.end) OVERLAPS (b.start, b.end);

The resulting file looks like this (purple lines):

Some borders are complete, others are interrupted (e.g. the border between Germany and the Netherlands) and many other borders are missing.

I think that some of the original polygons did not align perfectly and tried to resolve this using ST_Snap (see below). However this did not change the result:

    CREATE TABLE t.snap AS SELECT a.id, a."name", a."start", a."end",  ST_Snap(a.geom, b.geom, ST_Distance(a.geom,b.geom)*1.01) as geom FROM t.test As a, t.test As b WHERE ST_Intersects(a.geom, b.geom) AND a.id!=b.id AND (a."start",a."end") OVERLAPS (b."start", b."end");

Does someone have an how to fix this?

P.S. It seems that the "integrate" tool in ArcMap does a good job fixing this issue, but I was hoping to find a solution in PostGIS or R.


texts: Mapping election results using QGIS losing data?

texts:
I tried mapping a table that does not have a spatial geometric component (CSV file) into a shapefile in QGIS with a common ID column. But I noticed that most of the data which is on the CSV file does not show when I open the attribute table in the shapefile.

My CSV file has 728 rows and my shapefile has only 160 rows

After I join these two tables with a common ID column, only 160 rows are being shown on the shapefile attribute table.

Does anyone know how to combine all the data in the CSV to the shapefile?


texts: How to create two layers with different criteria from a raster layer

texts:
I have one raster layer I also turned into a slope layer. I need to create two rasters with different criteria ( first one needs specific slope and elevation, second needs just elevation) what tool would I use for this and how would I put it in the tool?
1 layer is elevation
I took that and made 1 slope layer
I need two rasters
Raster1: slope less than or equal to 5 and elevation less than 200
Raster2: elevation between 100 and 140


texts: Online service to transform coordinates from WGS84 to NZTM and NZMG?

texts:
Has anyone come across a free web service / REST API to transform coordinates?

I'm particularly interested in transforming from WGS84 to: New Zealand Geodetic Datum 2000 (NZGD2000 / NZTM) and New Zealand Geodetic Datum 1949 (NZGD1949 / NZMG). Would be a bonus if I could also get the grid version of coordinates as well.

Perhaps there's a JavaScript / C# library out there that I can use and host my own service?


texts: GRASS GIS - Keep raster below/above surface

texts:
I would like to clip raster where the surface intersect another surface and keep elevations that are below only.  I am interpolating a vector into a 2D raster (v.surf.rst) for properties at a given elevation (ex: 120 masl), and i would like to clip the raster when it is above the ground surface using a DEM.
A- I could subtract the DEM elevation from the raster elevation, and keep only positive values, this way only pixels which are below the DEM would be kept. Using R.Mapcalc.simple, you cannot select which field to apply the formula.


texts: Is there a way to use the old command line in Arcgis10?

texts:
I know that the command line has been replaced by the python window, but I was just wondering if there is any way to use command line within Arc 10?
 Thanks


texts: Using keyword in SelectByAttribute expression?

texts:
I've used SearchCursor to return the single value of a record (8 digit number)

for row in arcpy.da.SearchCursor(table, ['PID']):
if row[0] is not None:
        ID = row[0]
    break

Now I would like to use 'ID' in a query.

arcpy.SelectLayerByAttribute_management(Table2, "NEW_SELECTION",'ID' == ID)

This is not a valid expression and I'm not sure how to make the keyword ID usable.


texts: Convert CSV file format to Breakline CSV format file

texts:
I have a csv file with almost 36000 rows.The format of csv file is
x,y,Polygon number, Z

I would like to arrange this csv file in the below mentioned format. Let’s say in this example, Polygon number 2,3, 4, 5, points are in the 3rd column. I would like to take x,y,z of the each polygon vertices and write 3 before it starts that polygon vertices and END at the end of the polygon vertices. Here you will observe that polygon vertices 2 is with in number "3" and "END" in the output file.

Original CSV file

Export Breakline CSV format file


texts: How to import and patch multiple ASTER files in QGIS?

texts: I'm a new user and just began with GRASS, where I learned to import 30x30 resolution ASTER files and then patch them (the six I'm working with) into one map. On QGIS, it doesn't appear that there's a similar function and all six tiles are kinda staggered regarding the color alignment. THen, how do you save the patched file and change the color schemes? Lastly, how to change the lat-long to UTM? And all this, from the same GeoTIFF file?


texts: How to automate locally QGIS plugin building steps using exe or script?

texts: I have made a plugin name Area Calc on my local computer. The plugin work very well. I have followed  this. I don't want to release plugin. Is there any way to automate the whole above task through a script or exe which can follow all the steps required to build plugin on another computer. In short I would like to create a setup for my QGIS plugin.


texts:
I was working with LiDAR in Opengeo suite while following the Boundless tutorials (http://workshops.boundlessgeo.com/tutorial-lidar/)

In my recent post  titled 'LiDAR and Opengeo', I mentioned about the error with PDAL execution. I solved the issue by  editing the 'laz2pg.xml' file and added appropriate port, localhost and password in the xml file.

And now pdal commands can send request to read files. On the terminal, the white cursor blinks for a while and then stopped. At this point  the terminal no longer respond. Here is the error message:

rblzaman@ubuntu:~/Downloads/tileindex$ pdal pipeline laz2pg.xml
Requested to read 11418635 points
Requested to write 11418635 points
0Killed
rblzaman@ubuntu:~/Downloads/tileindex$

In the following day, in my same PC (Ubuntu installed on windows 7 via VM player),

it show the following error:

rblzaman@ubuntu:~/Downloads/tileindex$ pdal pipeline laz2pg.xml
Requested to read 11418635 points
Requested to write 11418635 points
0Caught exception: std::bad_alloc
rblzaman@ubuntu:~/Downloads/tileindex$

When I typed   'free -m' , it returns,
rblzaman@ubuntu:~$ free -m
             total       used       free     shared    buffers     cached
Mem:          1000        744        256          0         16        233
-/+ buffers/cache:        493        506
Swap:         1020        329        691
rblzaman@ubuntu:~$

Is this memory issue ? I have installed Ubuntu 12.04 on windows 7 using VM player. How to solve this error ?

Thanks


texts: What type of remote sensing methodology is commonly used to locate and track the movement of woodland caribou herds?

texts: I am most curious whether it is possible to use Landsat or other freely available imagery (but info on any other method would great as well) to locate and track the movement of woodland caribou herds over time.   I only have some experience doing this type of work with forest depletions and vegetation analysis but am curious whether similar techniques exist for locating and tracking animal herds.


texts: What is the absolute vertical accuracy and the relative vertical accuracy of SRTM 1-arc second?

texts: When I was looking for the accuracy values for the srtm dem, I found research papers trying to find out these values for its own selected area of study, some of them assume that the absolute value in Africa is around 5.6. I am new to studying remote sensing and GIS, so knowing where exactly to go for trusted information is a bit of challenge.