texts: Slope rasters - different results from different tools

texts:
I'm calculating Least Cost Path (LCP) within QGIS 2.18. The first step is to produce the slope raster from the DEM. Using r.slope in GRASS produces excellent results where I can go on and finally produce my LCP, whereas using the slope tool within QGIS or GDAL produces something quite different that does NOT allow me to produce my LCP. As a newbie I find this baffling. Can someone explain please? Surely these tools should at least produce something similar??

gdalinfo "E:/BU Volunteer/QGIS/218/gdal-slope.tif"

Driver: GTiff/GeoTIFF
Files: E:/BU Volunteer/QGIS/218/gdal-slope.tif
Size is 1824, 1416
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["unknown",
        SPHEROID["WGS84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (32.899999999999999,15.000000000000000)
Pixel Size = (0.008333333000000,-0.008333333000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  32.9000000,  15.0000000) ( 32d54' 0.00"E, 15d 0' 0.00"N)
Lower Left  (  32.9000000,   3.2000005) ( 32d54' 0.00"E,  3d12' 0.00"N)
Upper Right (  48.0999994,  15.0000000) ( 48d 6' 0.00"E, 15d 0' 0.00"N)
Lower Right (  48.0999994,   3.2000005) ( 48d 6' 0.00"E,  3d12' 0.00"N)
Center      (  40.4999997,   9.1000002) ( 40d30' 0.00"E,  9d 6' 0.00"N)
Band 1 Block=1824x1 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

gdalinfo "E:/BU Volunteer/QGIS/218/grass-slope.tif"

Driver: GTiff/GeoTIFF
Files: E:/BU Volunteer/QGIS/218/grass-slope.tif
Size is 1824, 1416
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (32.899999999999999,15.000000000000000)
Pixel Size = (0.008333332999970,-0.008333333000039)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.4.1 with GDAL 2.2.4
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  32.9000000,  15.0000000) ( 32d54' 0.00"E, 15d 0' 0.00"N)
Lower Left  (  32.9000000,   3.2000005) ( 32d54' 0.00"E,  3d12' 0.00"N)
Upper Right (  48.0999994,  15.0000000) ( 48d 6' 0.00"E, 15d 0' 0.00"N)
Lower Right (  48.0999994,   3.2000005) ( 48d 6' 0.00"E,  3d12' 0.00"N)
Center      (  40.4999997,   9.1000002) ( 40d30' 0.00"E,  9d 6' 0.00"N)
Band 1 Block=1824x1 Type=Float32, ColorInterp=Gray
  Description = slopeccfd2dd742084285a4245fec4de42bcc
  NoData Value=nan
  Metadata:
    COLOR_TABLE_RULES_COUNT=7
    COLOR_TABLE_RULE_RGB_0=0.000000e+000 2.000000e+000 255 255 255 255 255 0
    COLOR_TABLE_RULE_RGB_1=2.000000e+000 5.000000e+000 255 255 0 0 255 0
    COLOR_TABLE_RULE_RGB_2=5.000000e+000 1.000000e+001 0 255 0 0 255 255
    COLOR_TABLE_RULE_RGB_3=1.000000e+001 1.500000e+001 0 255 255 0 0 255
    COLOR_TABLE_RULE_RGB_4=1.500000e+001 3.000000e+001 0 0 255 255 0 255
    COLOR_TABLE_RULE_RGB_5=3.000000e+001 5.000000e+001 255 0 255 255 0 0
    COLOR_TABLE_RULE_RGB_6=5.000000e+001 9.000000e+001 255 0 0 0 0 0


texts: Merging three linestring tables into one new table view

texts:
I have three tables in a PostGIS database and I want to create a table view merging the features from the three tables into one, by using PostgreSQL.
The three tables are containing three different kinds of roads, hence I want the view to specify the type of road in a new column.

The three tables do all contain a geometry column (linestrings)

How can I do this?


texts: "Add Comment" functionality in Flex Viewer

texts: I want to add some kind of commenting functionality for users to leave notes on a map. I know that Esri deprecated the "add pop up" element in 2.3.1 because it didn't save. My current thinking is to direct users to create a circle feature and add their comment as an attribute using the "Edit" widget, however that seems needlessly complicated. Does anyone have a simpler solution for saving location-aware comments on the map?


texts: I processed MODIS MOD13Q1 NDVI times series (from 2000 to 2015) in Timesat. I specifically extracted start of season data and for each year (growth season) I was able to produce an image showing the SOS across my study area using the seas2img function. I was able to save these images as pictures but fail to understand how these images can be exported in a format that I will be able to use in ArcMAP in order to do further analyses. Please advise if this is at all possible and what process I need to follow to get it done.


texts: Move raster image to overlay it perfectly on another same resolution raster

texts:
This question may seem strange at first, but what I need to do is to move a raster image over another raster image and align them perfectly. Both rasters have the very same size and resolution - one is a Landsat band and another is a product of bands analysis, so this movement includes change of coordinates.

Why are they in different places? In Correcting Landsat images shifted by strangely high distance using ArcMap? I said I have rasters straight from Earth Explorer covering some area that are misplaced from 100  m even up to 40 km (don't ask me why, sensor callibration?).

I will explain it by listing what I do:

For every image create an RGB composition (for better terrain recognition)
Choose one image as a reference one (which seems to be correct according to my GPS data)
Using ArcMap's Shift tool I move each RGB image individually to it's correct place looking at the reference raster (RGB product only, bands are left untouched)
Now I would like to automatically apply movement of each band of each image to corresponding RGB image created from these bands (move them to their "correct place")
Afterwards I would use Snap raster option to align pixel misplacement if needed

Point 4. is where I'm stuck. I would prefer not to have to use Shift tool for each image, because selecting correct X and Y shift is problematic if you don't remember by how much you have already moved the corresponding RGB.


texts: Getting nearest Hospital information by using MapQuest API

texts:
How to get a list of hospital details (with address) near by customer address using MapQuest?

For example my address is : 1015 museum circle jacksonville FL

Using 5-box address format KVP Request, I got
Lat,Long
("lat":30.318548,"lng":-81.659093)

How to format the MapQuest or Nominatim URL request to get hospital list near by customer address?

Keyword to search : Hospital Or Clinic

Radius : 50 miles or 100 miles


texts:
I need to add Google Map with about 400 markers into WordPress site (into one page). Every marker needs to have info window with photo, description and link. It is really high number of markers, so I'm not sure what is the best way to implement that kind of map.

What do you suggest?


texts: QGIS 2.10 Can't load georeference plugin?

texts:
Apparently my problem was already present in older versions of Qgis. I don't know how to fix it, knowing I have a Linux system (Fedora 23).
Somethiing similar here :
Where to download QGIS 2.10?

So, I cannot found this plugin, neither Fetch Python Plugins command, but not it would fix anything anyway since georeference is within GDAL. I had installed  Qgis on this system (fedora 23) for first time (no older version).

Otherwise the procedure to download the plugin should be that one : http://www.qgistutorials.com/en/docs/georeferencing_basics.html
but I do not find georeference into the list of available plugins. The repository I use is the default one : http://plugins.qgis.org/plugins/plugins.xml


texts: Calculating area gets a different result

texts:
I've calculated the area of polygons with QGIS by using ctr-i - geometry - $area. I've set it on real. The results seem fine, but I've noticed they differ from the area given when I use the i-function.

For example: polygon 4 gives an area of 7,662 km2 when using the i-function but when I Calculate it, it gives an area of 7682816.27364349.

Does anybody know what I'm doing wrong?

Layer CRS is Amersfoort / RD new EPSG:28992


texts:
I'm using Openlayers 4.6 API and the overviewMap function is not working.
I'm using Bingmap areial as basemap.

  this.map = new ol.Map
  ({
      controls: [new ol.control.OverviewMap()],  //, attr
      interactions: ol.interaction.defaults
      ({
          doubleClickZoom: false
      }),
      layers: [Bing_layer],
      target: 'map',
      view: view
  });


texts: Getting Null Values After "toreal" Command

texts:
I  have joined a simple CSV file to my project and am attempting to change a string value to a real value.  In editing mode, I bring up the field calculator and create a new field from the first field.  I try to change the string to real by doing the command "toreal(nameofmyfield)".  When I create the new field though, I get NULL for every value, even though I see the normal numbers in the column I am trying to change.

Any idea what's going on?


texts:
I'm attempting to export the current view to a TIFF (I attempted JPEG as well) image in ArcGlobe, but it seems most paths to this are broken in ArcGlobe even though the first two of these work in ArcMap. TIFFs can be exported in ArcGlobe within the application menu so some way must work. Is there another way?

This one gives me a NotImplementedException

    void CreateTIFFFromActiveView(IActiveView activeView, string pathFileName)
    {
        //get pixel boundary of the active view
        Rect rect = default(Rect);
        GetWindowRect(activeView.ScreenDisplay.hWnd, ref rect);
        var envelope = new Envelope() as IEnvelope;
        envelope.PutCoords(rect.Left, rect.Top, rect.Right, rect.Bottom);
        var exportRect = new tagRECT()
        {
            top = rect.Top,
            bottom = rect.Bottom,
            right = rect.Right,
            left = rect.Left
        };

        IExport export = new ExportTIFF() as IExport;
        export.ExportFileName = pathFileName;
        export.Resolution = 600;
        export.PixelBounds = envelope;
        int hDC = export.StartExporting();

        //NotImplementedException HERE
        activeView.Output(hDC, 600, ref exportRect, activeView.Extent, null);

        // Finish writing the export file and cleanup any intermediate files
        export.FinishExporting();
        export.Cleanup();
    }

This one causes an AccessViolationException

    void exportTif(string outputPath)
    {
        var activeView = ArcGlobe.Document.Scene as IActiveView;

        Rect rect = default(Rect);
        GetWindowRect(activeView.ScreenDisplay.hWnd, ref rect);
        var envelope = new Envelope() as IEnvelope;
        envelope.PutCoords(rect.Left, rect.Top, rect.Right, rect.Bottom);

        var docPrintExport = new PrintAndExport();

        var tifExport = new ExportTIFF() as IExport;
        tifExport.ExportFileName = outputPath;
        tifExport.PixelBounds = envelope;
        tifExport.Resolution = 600;

        //AccessViolationException
        docPrintExport.Export(activeView, tifExport, 600, false, null);
      }

This one just leaves me with a blank image

    void drawBmp()
    {
        _bmp = new Bitmap(2000, 1000);
        _gfx = Graphics.FromImage(_bmp);
        _hdc = (int)_gfx.GetHdc();

        var globeDoc = ArcGlobe.Application.Document as IGMxDocument;
        var globe = globeDoc.Scene as IGlobe2;
        var activeView = globe as IActiveView;
        activeView.Draw(_hdc, null);

        //just wait 30 seconds for drawing and whatever else
        var timer = new System.Timers.Timer(60000);
        timer.Elapsed += (o, args) =>
        {
            _context.Post(x =>
            {
                _bmp.Save(outputPath);
                _gfx.ReleaseHdc((IntPtr)_hdc);
                _gfx.Dispose();
                _bmp.Dispose();
            }, null);
        };
        timer.AutoReset = false;
        timer.Start();
    }

This one works, but is limited to the current resolution of the active view which is smaller than I need

    void snapShot(string outputPath)
    {
        var globeDoc = ArcGlobe.Application.Document as IGMxDocument;
        var globe = globeDoc.Scene as IGlobe2;
        var activeView = globe as IActiveView;

        Rect rect = default(Rect);
        GetWindowRect(activeView.ScreenDisplay.hWnd, ref rect);

        globe.GlobeDisplay.ActiveViewer.GetSnapshot(rect.Right - rect.Left, rect.Bottom - rect.Top, esri3DOutputImageType.BMP, outputPath);
    }


texts: Create an X,Y latitude and longitude grid from an raster

texts:
I have a satellite image (LandSat 5), I want create a raster with X,Y coordinates, this is possible in Model Builder of ArcGis 10.0? or other option. I know that the option $$YMAP ($$XMAP) not work more in ArcGis 10.0

I want make the same. But my input is a raster (landsat image, UTM projection), I need obtain (output) a grid  of latitude and a grid of longitude (two grids: latitude and longitude). First, this  code is valid, and if, is possible work with ModelBuilder. I'm novel in arcgis and python, I working with evapotranspiration and I need this grids for calculate it.
This is the code of video (www.youtube.com/watch?v=_-6apNB18gk):


texts: What is the procedure for reflectance calculation from Sentinel 2 1C products downloaded from ESA Scientific Hub?


texts: I was wondering if anyone could help me out with a website which I can use to get geometries for a area. What I'm looking for is something like google maps where I can drag a cursor over a area of the map and get a bounding box for it (as geometry not geography). If anyone knows of a site like this or one where i can find out bounding box's it would be greatly appreciated.


texts:
I was happily using html popups with QGIS 2.8.3 for various things including a showing a photo. However, I've updated to 2.18 to take advantage of other new functionalities, and now find that my html popup will no longer show the photo the way it used to.

This is the html I'm trying to use:

The image name is in field PHOTO_1 and the image is in a subfolder named Camera (subfolder on same level as qgis project).

With version 2.8.3, the html provided showed the photo. With version 2.18, I only have a blank popup.

I've now tried it two additional ways without success:

and


texts: Remove Interactions and Clear Vectors Openlayers3?

texts:
I have a simple mapping application that allows a user to select a draw type eg polygon or point and then draw shapes on the map.

I have to use buttons instead of a dropdown list to select the draw type, my issue comes when I try to clear the previous draw type, eg the user clicks polygon > draws shape > changes their mind > clicks point > map sucessfully clears previous interactions > but when they go to draw the shape now draws both a polygon and a point.

Below is a snippet of code for the draw point event trigered by a button press

        document.getElementById("drawPoint").addEventListener("click", function drawPoint() {
      var geometryFunction, maxPoints;
      draw = new ol.interaction.Draw({
        source: source,
        type: ("Point"),
        geometryFunction: geometryFunction,
        maxPoints: maxPoints
      });
      vector.getSource().clear();
      map.addInteraction(draw);
      source.on('addfeature', function(evt){
          var feature = evt.feature;
          var coords = feature.getGeometry().getCoordinates();
      });

    });

the polygon button is the exact same apart from shape etc.
I am trying to use the following function to clear the map and the interactions but seem to be missing something.

        document.getElementById("myForm").addEventListener("change", changeFunc);

    function changeFunc() {
        map.removeInteraction(draw);
        vector.getSource().clear();
        addInteraction();
    }


texts: How Can I Install "QGIS Server" in a Windows 7 PC?

texts:
I am trying to install the new "QGIS Server" feature that comes with QGIS 1.6 on my Windows PC.  The directions for this at the QGIS website only show instructions for Mac OS, and those in other GIS blogs show instructions for installing the feature on Linux machines.

I'm wondering if instructions for installation on Windows are available anywhere, or whether anybody can provide me with the steps for installing the feature on the Apache web server on my PC.


texts: Geonetwork 2.10: creating metadata for documents (reports, etc.)

texts:
I am using geonetwork 2.10, and I would need to create meta data for documents such as reports.

Which template am I supposed to use?
The only template I found which looks close enough to me is "template for map - Use this template to describe a static map (eg. PDF or image) or an interactive map" from ISO 19139.

Can you suggest me another template?


texts: Re-projection issues with Azimuthal Equidistant

texts:
When reprojecting country-border polygons from longlat:

+proj=longlat +datum=WGS84 +no_defs

to Azimuthal Equidistant:

+proj=aeqd +lat_0=26 +lon_0=126 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

I occasionally find a polygon that is turned inside out. The settings above show an example of this. The first image shows the polygon for Argentina plotted 'as-is', and the second image shows the same data that has had long lines densified and then re-projected. The scale is adjusted for clarity. The 3rd image shows All Countries, and their polygons would be filled correctly except that the inside-out polygon for Argentina floods them all. How can I determine when a polygon will display 'inside out', and how can I then make it display correctly? I am using DotSpatial to do the re-projection, but other methods produce the same results. A grid has been overlaid for clarity.


texts: Is it possible to programmatically add calculated fields?

texts:
In the Memory provider described in PyQGIS Cookbook or in question How to create a new empty vector layer programmatically? I can see how to create a vector layer programmatically using standard fields (String, Int, Double).

with PyQGIS/QGIS 1.8.0, is it possible to programmatically add calculated fields like the ones defined using the GUI field calculator? After some reading about QgsExpression I'm thinking in something like this:

exp = QgsExpression('Column * 3')
pr.addAttributes( [ QgsField( "fieldName", QVariant.Expression, exp ) ] )


texts: Choosing projection to use for Least-Cost Path analysis at continental scale?

texts: I am trying to find least cost paths from the southern USA to the Great Lakes region (for example). My resistance layers are currently in an equal area projection but I'm wondering if I should convert to an equidistant projection? My understanding is that part of what factors into a least-cost path analysis is distance (all else being equal, the shortest route from A to B is the least-costly) and so to me it would make sense to have the raster layers in an equidistant projection...but maybe I'm missing something? If equidistance does make sense, what is the most appropriate for all of North America?


texts: QGIS expressions function not loading when opening composer

texts:
I created a function to display the coordinate system of a map and transformation parameters etc. When I open the composer, I have to go to functions and press on load button to re-load my function to make it work. Otherwise, I get an error as follows:

Eval Error:
global name 'qgis' is not defined

Code is as follows:

"""
Define new functions using @qgsfunction. feature and parent must always be the
last args. Use args=-1 to pass a list of values as arguments
"""
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
import re

@qgsfunction(args='auto', group='Custom')
def test5(feature, parent):
    canvas = qgis.utils.iface.mapCanvas()
    crsname = canvas.mapRenderer().destinationCrs().description()
    return crsname

I am using 2.18.14 of QGIS. As soon as I hit the load button, the code starts working but I want to make sure it auto-loads on startup.


texts: How to show elevations within a 5% slope or less in ArcMap 10.4?

texts: I have a vector contour shapefile that shows elevation lines (i.e. USGS topo contours). I want to highlight areas that have a slope of 5% or less. I was wondering if someone could explain how to do this using ArcMap 10.4? I thought building a Query could work, but I'm not sure where to start with that.


texts: Why pgr_dijkstra function is faster with directed graph?

texts:
I'm doing thousands of tests to get the average route time calculated with the pgr_dijkstra function, taking random nodes.
In all tests the average time to generate the routes of the directed graphs (vehicles) is about 100 milliseconds faster than for undirected graphs (pedestrians). It seems like a very short time, but generating thousands of routes the total time less turns into minutes.
Why is that? Is it because in the directed graphs the reverse_cost is ignored, and that is why the function does not spend its time analyzing this field? And in the undirected graphs is the cost and reverve_cost field analyzed?
Function example with undirected graph:
SELECT * FROM pgr_dijkstra(
'SELECT gid AS id,
     source,
     target,
     length AS cost
    FROM ways',
9411, 3986,
directed := false);


texts: Splitting an administrative region into given number of parts

texts:
I have the following plot.

I want to divide the administrative region into given number (say 20) equal parts.

library(maps)
data(world.cities)

Pakistan <- data.frame(map("world", "Pakistan", plot = FALSE)[c("x","y")])

library(ggplot2)

p <- ggplot(Pakistan, aes(x = x, y = y)) +
     geom_path(colour = 'green', linestyle = 2) +
     coord_map() +
     theme_bw() +
      labs(x=" ", y=" ") +
  theme(
         panel.grid.minor = element_blank()
      , panel.grid.major=element_blank()
      , axis.ticks = element_blank()
      , axis.text.x = element_blank()
      , axis.text.y = element_blank()
      , panel.border = element_blank()
      )

print(p)

Edited

I found this one

makeVchopper <- function(pol){
    bb = bbox(pol)
    delta = (bb[2,2] - bb[2,1])/10
    xmin = bb[1,1]-delta
    ymin = bb[2,1]-delta
    ymax = bb[2,2]+delta

    choppoly = function(xmax){
        readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
                        xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin, xmin,ymin))
    }
    choppoly
}

slicer <- function(pol, xmin, xmax){
    bb = bbox(pol)
    delta = (bb[2,2] - bb[2,1])/10
    ymax = bb[2,2] + delta
    ymin = bb[2,1] - delta
    r = readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
        xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin, xmin,ymin))
    gIntersection(pol,r)
}

chop_thirds <- function(pol, fractions=c(1/3, 2/3)){
    chopper = makeVchopper(pol)
    bb = bbox(pol)
    xmin = bb[1,1]
    xmax = bb[1,2]

    totalArea = gArea(pol)

    chopped_area = function(x){
        gArea(gIntersection(chopper(x),pol))
    }

    edges = lapply(fractions, function(fraction){
        target = totalArea * fraction
        target_function = function(x){
            chopped_area(x) - target
        }
        uniroot(target_function, lower=xmin, upper=xmax)$root
    })

    xdelta = (xmax-xmin)/10
    chops = matrix(c(xmin-xdelta, rep(edges,rep(2,length(edges))),
xmax+xdelta), ncol=2, byrow=TRUE)
    apply(chops, 1, function(edges){
        slicer(pol, edges[1], edges[2])
    })

}

# Usage:

library(rgeos)
library(sp)

# sample data
pol <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20), (-150 -20, -100 -10, -110 20, -150 -20))")
plot(pol)

# now split

parts = chop_thirds(pol)
plot(pol)
plot(parts[[1]], add=TRUE, col=1)
plot(parts[[2]], add=TRUE, col=2)
plot(parts[[3]], add=TRUE, col=3)

gArea(parts[[1]])
gArea(parts[[2]])
gArea(parts[[3]])


texts: Associate output of QGIS GRASS distance/network functions with 'from' points rather than 'to' lines

texts:
I am calculating the distance from remote locations, which I have as point vector data, to a road network, a line vector. I'm using QGIS 2.18.2.

As per this image, I've had success with using the integrated GRASS tools v.distance and v.net.connect to generate the lines to the road network (though with slightly different results between the two functions, hence the green and red). However, the attribute file of the new line layers generated take their attribute data from the roads they're connecting to, rather than the locations they originate from.

QGIS can give me the distance represented by the generated lines, but I need to associate those distances with the locations.

Is there a different function, or a subsequent one I should be using? Ultimately I want to end up with the distance to the road as an additional attribute of the original points.


texts: New Geocoded Points Layer Misaligned With Polygon Layer

texts:
I learned on ArcGIS, and I'm familiarizing myself with the use of QGIS.  So far I love it, but the process of re-projecting and transforming datums seems less intuitive than in Arc.

I recently geocoded 63 locations, pulled in the csv and transformed into a points shapefile. I then set the projection for that points layer as the same for my block groups multi-polygon (EPSG 102121- NAD_1983_Michigan_GeoRef_Feet_US).

The 2 layers are still misaligned, even if I save each layer with the specified projection and try to reopen in a blank project.  Can anyone help?


texts:
I have a method that does a bunch of geoprocessing on some layers for my ArcMap extension. It creates a feature-to-point layer for one FeatureClass and another feature-to-point layer for another FeatureClass. Then one of the feature-to-point layer's is used to as input to create a Buffer layer. This all works without any errors until I try to run the method in a separate thread. When I run it in a separate thread, both feature-to-point tools run without errors but it fails when it tries to create the buffer.

I borrowed error handling for my geoprocessing from Avoiding fails from ArcObjects geoprocessing with .NET?. The result shows the following.

Executing: Buffer GPL0 C:\Users\me\Documents\geoprocess.gdb\buffer "250 Feet" FULL ROUND NONE # esriGPMessageTypeProcessDefinition
Start Time: Tue Mar 06 18:55:10 2012 esriGPMessageTypeProcessStart
Failed to execute. Parameters are not valid. esriGPMessageTypeError
ERROR 000732: Input Features: Dataset GPL0 does not exist or is not supported

It seems when I try to run this in a separate thread, the name of the input FeatureClass to the buffer tool gets lost and changed to "GPL0" which doesn't exist. I've read that there can be issues when passing objects between threads but all the geoprocessing is happening in the same thread (as far as I understand). I pass an IMap variable and some other objects to the method. The IMap object seems to be passed okay since I can get its layers and the other objects are not used for the geoprocessing.

In my form:

backgroundWorker.DoWork += (obj, eArgs) => BgWorkerDoWork(map, ...)
bgwCalibrate.RunWorkerAsync();

private void BgWorkerDoWork(IMap map, ...)
{
  var thread = new Thread(() => MyClass.MyGeoProcessMethod(map, ...));
  thread.SetApartmentState(ApartmentState.STA);
  thread.Start();
  thread.Join();
}

In my class MyClass:

private void MyGeoProcessMethod(IMap map, ...)
{
  // do stuff...
  // then create both feature-to-point files
  IFeatureClass inFeatureClass1 = map.Layer[0] as IFeatureClass;
  IFeatureClass inFeatureClass2 = map.Layer[1] as IFeatureClass;
  // no errors:
  IFeatureClass polygonsToPoints1 = CreateFeatureToPoint(inFeatureClass1);
  IFeatureClass polygonsToPoints2 = CreateFeatureToPoint(inFeatureClass2);
  // Create buffer (this is the part that fails)
  IFeatureClass bufferFeatureClass = CreateBuffer(polygonsToPoints2);
}

Both feature-to-point layers get created successfully. The geoprocessing events say so and I can see the results on the map. Any suggestions or ideas on what's going on with the buffer part?


texts: Save points and lines within the same shapefile

texts: Is it possible to save points and multiline strings within the same shapefile? If so, how?


texts:
I try to open a Lidar file in GRASS GIS 7.8.5 and 7.9 in a Macintosh computer. However, when I try to open the tool  r.in.lidar, it appears:
No such file or directory: 'r.in.lidar'.


texts:
Sometimes when trying to construct an sf object from an sfc object I get the following error:

“Error in MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) :polygons not (all) closed”.

I am sure that my polygons are closed and I strongly assume that the floating error problematic causes this error. I thought a solution might be to round all coordinates of the sfc object to 1m (assuming that I am using a cartesian coordinate system). I have, however, not succeeded. Any thoughts?

library(sf)
p1 <- rbind(c(1000.001,1000.002), c(1000.2345,2000.23452), c(2000.23345,2000.2342), c(1000.001,1000.002))
p2 <- rbind(c(500.234,500.3234), c(500.234,100.235), c(1000.00001,1000.234), c(500.234,500.3234))
pol <-st_polygon(list(p1, p1))
pol_sfc= st_sfc(pol)
pol_sf <- st_sf(pol_sfc)
st_coordinates(pol_sf)


texts: Problems with simple Leaflet + Mapbox TitleLayer?

texts:
I'm completely new to GIS, and thought I would start off with the leaflet  quick tutorial. However, I'm having problems with adding the base tile layer with Mapbox. This is my code so far:

    <!DOCTYPE html>

    <title>Agricultural Trade Patterns</title>
    <meta charset='utf-8' />
    <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no" />
    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />

I tried to follow the tutorial as much as I could, but whenever I test it I'm getting the error:

TypeError: L.titleLayer is not a function

And I can't figure out why.

Am I passing the wrong arguments to tileLayer?

I know my URL template doesn't have the exact format that was on the leaflet tutorial, but I couldn't find a png url like that anywhere?

Is that the problem?


texts:
Is there an issue with the SAGA-tool "Mosaicking rasters"?
I tried it with different rasters. In Sage itself it works fine.

  Algorithmus Mosaick raster layers startet...
  io_gdal 0 -TRANSFORM -INTERPOL 0 -GRIDS "C:\Temp\processing\04593242dc054133b9fbcae916bae1f8\1tif.sgrd" -FILES "C:\Temp\test\2.tif,C:\Temp\test\1.tif"
  grid_tools "Mosaicking" -TARGET 0 -GRIDS "C:\Temp\processing\04593242dc054133b9fbcae916bae1f8\1tif.sgrd" -TYPE 7 -INTERPOL 0 -OVERLAP 1 -BLEND_DIST 10 -MATCH 0 -USER_FIT 0 -USER_XMIN 331858.3 -USER_XMAX 335568.9 -USER_YMIN 5921607.7 -USER_YMAX 5924629.6 -USER_GRID "C:\Temp\processing\0aab6cedcd08418eab5e1ee4eca47901\USERGRID.sdat"
  C:\PROGRA~1\QGISPI~1\bin>set SAGA=C:\Program Files\QGIS Pisa\apps\saga
  C:\PROGRA~1\QGISPI~1\bin>set SAGA_MLB=C:\Program Files\QGIS Pisa\apps\saga\modules
  C:\PROGRA~1\QGISPI~1\bin>PATH=PATH;C:\Program Files\QGIS Pisa\apps\saga;C:\Program Files\QGIS Pisa\apps\saga\modules
  C:\PROGRA~1\QGISPI~1\bin>saga_cmd io_gdal 0 -TRANSFORM -INTERPOL 0 -GRIDS "C:\Temp\processing\04593242dc054133b9fbcae916bae1f8\1tif.sgrd" -FILES "C:\Temp\test\2.tif,C:\Temp\test\1.tif"
  ERROR 1: Can't load requested DLL: C:\Program Files\QGIS Pisa\apps\saga\dll\gdal_MrSID.dll
  127: Die angegebene Prozedur wurde nicht gefunden.
  ERROR 1: Can't load requested DLL: C:\Program Files\QGIS Pisa\apps\saga\dll\gdal_MrSID.dll
  127: Die angegebene Prozedur wurde nicht gefunden.
  ERROR 4: C:\Temp\test\2.tif,C:\Temp\test\1.tif' does not exist in the file system,
  and is not recognised as a supported dataset name.
  library path: C:\Program Files\QGIS Pisa\apps\saga\modules\io_gdal.dll
  library name: Import/Export - GDAL/OGR
  tool name : GDAL: Import Raster
  author : O.Conrad (c) 2007 A.Ringeler)
  Parameters
  Grids: No objects
  Files: "C:\Temp\test\2.tif,C:\Temp\test\1.tif"
  Transformation: yes
  Interpolation: Nearest Neighbor
  loading: C:\Temp\test\2.tif,C:\Temp\test\1.tif
  failed: could not find a suitable import driver
  Error: Dataset creation failed:C:\Temp\test\2.tif,C:\Temp\test\1.tif' does not exist in the file system,
  and is not recognised as a supported dataset name.
  Error: executing tool [GDAL: Import Raster]
  C:\PROGRA~1\QGISPI~1\bin>saga_cmd grid_tools "Mosaicking" -TARGET 0 -GRIDS "C:\Temp\processing\04593242dc054133b9fbcae916bae1f8\1tif.sgrd" -TYPE 7 -INTERPOL 0 -OVERLAP 1 -BLEND_DIST 10 -MATCH 0 -USER_FIT 0 -USER_XMIN 331858.3 -USER_XMAX 335568.9 -USER_YMIN 5921607.7 -USER_YMAX 5924629.6 -USER_GRID "C:\Temp\processing\0aab6cedcd08418eab5e1ee4eca47901\USERGRID.sdat"
  ERROR 1: Can't load requested DLL: C:\Program Files\QGIS Pisa\apps\saga\dll\gdal_MrSID.dll
  127: Die angegebene Prozedur wurde nicht gefunden.
  ERROR 1: Can't load requested DLL: C:\Program Files\QGIS Pisa\apps\saga\dll\gdal_MrSID.dll
  127: Die angegebene Prozedur wurde nicht gefunden.
  library path: C:\Program Files\QGIS Pisa\apps\saga\modules\grid_tools.dll
  library name: Grid - Tools
  tool name : Mosaicking
  author : O.Conrad (c) 2003-12
  Load grid: C:\Temp\processing\04593242dc054133b9fbcae916bae1f8\1tif.sgrd...
  failed
  Mosaicking: could not initialize data objects
  Error: Grid file could not be opened.
  Error: executing tool [Mosaicking]
  C:\PROGRA~1\QGISPI~1\bin>exit
  Converting outputs
  Lade Ergebnis Layer


texts:
Is it possible to resample the raster by using GDAL C# binaries, in a way that every other cell in a row and column is removed?

Basically a replication of this GDAL command:

gdal_translate -outsize 50% 50% C:/input.tif C:/output.tif

Looks like I need to create a new dataset first:

outputDataset.GetDriver().Create("C:/output.tif", inputDataset.RasterXSize/2, inputDataset.RasterYSize/2, channelCount, DataType.GDT_Byte, null);

Then read the elevation values from the initial (non-resampled) dataset by applying the buffer size to be half of those the inputDataset had (inputDataset.RasterXSize/2 , inputDataset.RasterYSize/2):

inputDataset.ReadRaster(0, 0, inputDataset.RasterXSize, inputDataset.RasterYSize, buffers, inputDataset.RasterXSize/2, inputDataset.RasterYSize/2, inputDataset.RasterCount, [inputDataset.RasterCount], 0, 0, 0);

And finally write the elevation values to a new .tif file:

finalDataset.WriteRaster(0, 0, inputDataset.RasterXSize, inputDataset.RasterYSize, buffers, inputDataset.RasterXSize/2, inputDataset.RasterYSize/2, inputDataset.RasterCount, [inputDataset.RasterCount], 0, 0, 0);

But the final .tif file I am getting does not seem to be correct.

Can somebody provide a C# example please?


texts:
Is there anyway to determine projection of a shapefile which doesn't have .prj file?

I am trying to re-project this shapefile to another projection, I know that it is Lat long data. I could manage to load this in GIS tools like this or QGIS

Based on DBF data, it's Auckland data. According to this link the EPSG code of New Zealand is 27120 as the picture indicates:

Also I tried to display a point of this shape in Google Maps and it worked:

But still not sure what is the projection of this shapefile? I don't think it's WGS84 because I wasn't successful to reproject it from WGS84 to EPSG:27120 which is New Zealand grid system.

If it's not WGS84 how does Google display the point successfully?

The shapefile is available here, it is small (5KB)


texts:
ITRF 96 1998.0 epoch

    X          Y           Z         Vx(m/y)   Vy(m/y)  Vz(m/y)
4059623.25 2764405.93 4056720.43    -0.0079   -0.0007  0.0132

Transformation parameters from ITRF1996 to ITRF2014.(epoch:2010)

Solut.  Tx      Ty      Tz     D       Rx        Ry        Rz
       (mm)    (mm)    (mm)  (ppb)   (,001)    (,001)    (,001)
Rates   Tx      Ty      Tz     D       Rx        Ry        Rz
       (mm/y) (mm/y) (mm/y) (ppb/y) (0,001/y) (0,001/y) (0,001/y)
ITRF96 -7.4     0.5    62.8  -3.80      0         0       -0.26
rates  -0.1     0.5    3.3   -0.12      0         0       -0.02

GRS 80 ellipsoid a=6378137.00 b=6356752.314m 1/f=298.257222101

Question is 2017,8911704 epoch ITRF2014 XYZ coordinates

For now I know how to convert epoch but I don't know how to convert datum(ITRF to ITRF). This is my quiz question. Is there anyone who can help me?


texts: How to convert a polygon to a line? essentially reverse buffering

texts:
So basically I've got this data that some other company sent me and instead of them sending it in line format like all the other companies, whatever troll works there just decided to send a 25meter buffer of the network

Is there a way to obtain the reverse of the buffer? Apparently ST_ApproximateMedialAxis would give me what I need but unfortunately amazon does not accept that extension. I'm also using QGIS


texts:
What is causing this strange behaviour in my map? As closer to the poles, the labels are distorted.

Local OpenStreetMap data in PostgreSQL
Geoserver
OpenLayers
Projection of data: 4326
Projection of map (OpenLayers): 4326

I don't know what more I can tell. Please ask me for more information.

Edit to add OpenLayers code asked by @Iant

        landLayer = new ol.layer.Tile({
            source: new ol.source.TileWMS({
                url: MCLM.Map.geoserverUrl,
                isBaseLayer : true,
                projection: ol.proj.get('EPSG:4326'),
                params: {
                    'LAYERS': MCLM.Map.baseLayerName,
                    'FORMAT': 'image/png8',
                    'tiled': true,
                    'antialias' : 'full',
                    'VERSION': '1.3.0',

                }
            })
        });


texts:
How can I find how many duplicates I have in my dataset (areas, lines, points) using ArcMap?

With topology I can find overlapping objects, but I want to find duplicates.

Are there any tools from toolbox?


texts:
I need to plot an inset on my raster map showing a zoomed-in part of the map (according to specificed coordinates). Something like the image below (credit: David Wachsmuth, twitter).

I know there are a lot of resources on how to do it in ggplot2 and shapefiles, but I've found nothing about how to do it with a raster map. A reproducible example of my raster is shown below:

 kiwi <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
 kiwi[] <- sample(2:35, ncell(kiwi), replace = T)

And I would like to create an inset of the following coordinates:

clipExt <- extent(-120, 10, 70, 90)
f<-crop(df, clipExt)


texts:
When I work with ArcGIS and georeference a picture, my picture (Tif) remains the same but new files are being created:

.tfwx .aux (XML) .ovr

The data of georeferenced picture are being stored within these files.

When I work wit Qgis, i get a modified picture and only the .ovr file.

I was wondering if there is a way to georeference in QGIS as it is being done in ArcGIS.


texts: Extract from raster by tightest rectangle mask without Nodata

texts:
I'm using ArcMap 10.2.2 and I have a raster that has Nodata around the edges.  How can I extract the tightest rectangle from this raster that does not include Nodata?  The screenshot may help illustrate my question.

This is my raster and Nodata is indicated as yellow.  What tool/script can I run/write that will extract the raster based on the red bounding rectangle?  The purpose of this is to NOT use manual/visual inspection to get the red bounding rectangle extent values.  I intend to use whatever method that works to include in a larger Python script.


texts: Show pie-chart values & legend in info-window for AGS_JSAPI 3.4/3.5

texts:
Currently I am going through this example and trying to insert values in pie chart (I have added that screenshot and want to show values & legend in the info window)

I am using ArcGIS Suit 10.1 and ArcGIS JSPAI version 3.5

What I have tried:

-Tried to insert values in media-info but it is not representing in pie-chart.

-Tried this dojo-example link

this is what I am trying to do in info-window :

Any help regarding this will be great

Thanks in advance : )


texts: MMQGIS: How to make grid layer in meters not degrees?

texts:
I've read several of the posts on here pertaining to MMQGIS Create Grid Layer option, but I have not found information on how to make the grid based on meters veres degrees. I know how to manipulate the CRS for shp layers to adapt to meters, but I'm not quite sure on how to do it within this context.

Specifically, I'd like to make a hex grid with roughly a radius of 804.672 meters for each hex. The grid would have to cover an area that is about a 48280.3 meter radius.

Am I doomed to inaccurate degree measurements of the Hex Grid, or is there a solution to this problem. I have a basic proficiency with QGIS. Thanks!


texts:
I no longer see the plugins, the processing toolbox and all the scripts in QGIS 3.10 LTR. They're not even in the "Processing" Menu. I tried looking to activate it from the menu bar but it's not the only problem, QGIS all in all is not connecting to the scripts.

I think I know that the problem is with python, I just need some help reversing it.

This is what I did:

I was experimenting with GDAL in th OSGeo4W shell
I needed to execute gdal_calc so following the tutorial I needed to set the path to the python app by exectuting the following code: py3_env
I executed py3_env many times after so it resulted in concatinating the paths as follows

When I launched QGIS 3.10 this is what I have

The not even a "Processing" Menu.

I imagine this requires that I reassign the path to the scripts ?

How can I fix this as I don't want to reinstall QGIS again.


texts: Generating MultiPoint column from two Point columns without changing number of rows using shapely

texts:
from shapely.geometry import MultiPoint, Point
a=Point(1, 1.5)
b=Point(2, 1.5)
c=MultiPoint([a,b])

I am trying to create a MultiPoint column from two Point columns. Using the example code above, I can obtain a multipoint c but the data length is doubled, which I cannot use for an entire dataset. I tried different brackets, none has worked.
My dataset has 42 rows and contains columns of centroid_x (Point(a,b)) and centroid_y (Point(a,b)), which are original and destination locations. When I creating a MultiPoint column that combining (centroid_x and centroid_y), it comes out with 84 rows. I tried list(zip(tuple(a), tuple(b))). It solved the row number issues, but couldn't work for MultiPoint. PS: the method I used works fine with LineString.


texts:
I am attempting to create a web map in OpenLayers by pulling in layers from GeoServer which is running on an Amazon EC2 instance.  The prepackaged layers from GeoServer seem to show up fine when I substitute their URL in the map code. However I have a PostGIS layer that I am wanting to use as the basis for the map which is not showing up no matter what I try.

I have read several threads that seem to indicate it's a problem with mismatched projections. But I'm not sure where or how to "reproject" the layer or if that is my problem. It does shows up in the GeoServer OpenLayers preview (as a single tile, but not "tiled" for some reason). The declared SRS when I originally published the layer in GeoServer is EPSG:404000. I'm not sure why it is that. Changing it doesn't seem to help.

<!DOCTYPE html>

    <title>Tiled WMS</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.4/css/ol.css" type="text/css">


texts: How to modify parameter file in MRTSwath for bulk data processing?

texts: I need to create a parameter file for bulk data processing in MRTSwath Tool. I have the basic script written but there is issue in modifying it so that it automatically picks the input files from a folder and processes it according to the given parameter.


texts:
Once again I encounter something very strange with ArcGIS. I perform a buffer analysis on one side of a polyline. There seems to be a serious accuracy problem and I can't figure out what the cause is. So I hope to find some expertise out there...What I tried without success:

change distance settings from 'Field' to 'Linear unit'
change distance settings from 'Field' to 'Linear unit'
change end type settings from 'ROUND' to 'FLAT'
run check geometry tool on the line feature
I made a copy of the line and shorted it to area of interest which had some effect but didn't resolve the problem entirely. Especially around the sharp angle the buffer still takes a shortcut!

The buffer is lurching up to 4 meters on the wrong side! See the picture I attached to get an better impression of the problem. Does anyone have an idea what might me causing this behavior? I'm using ArcGIS 9.3.1. with WinXP.

UPDATE

See also my comment below.

I performed the described series of actions (divide, intersect,..) to increase the accuracy level of the source data, placing a vertex every 3 meters along the line. There still remains inaccuracy, again the buffering around the rough turn especially is giving me trouble. See second pic.


texts:
I have produced a loop which step-wise applies a mask with increasing threshold to an image.  Using the Map.addLayer() feature, I wish to view all outputs however I do not seem to be able to name them different things within the loop:

for (var i = 0.85; i < 0.96; i += 0.01){

    //sets the threshold
    var threshold = ee.Number(i);

    //two input images from the minumum and the error
    var min = ee.Image(image4);
    var error = ee.Image(image3);

    //applies the threshold to the minimum image
    var percent_mask =
    min.select('b1').lte(error.select('b1').multiply(threshold));
    var masked_percent = min.updateMask(percent_mask);

    Map.addLayer(masked_percent, {}, 'name');

}

I wish to replace 'name' with something unique, preferably something to do with the threshold applied for each output.
Full code with assets can be found here:
https://code.earthengine.google.com/6e4c1c3c4691e29f7149176447076197


texts: Updating Overpass-API database with local changes file

texts:
I have a local OSM server and its data is stored in openstreetmap PostgreSQL database. Also I have Overpass-API installed for the local OSM server (using this tutorial).

The openstreetmap database is sometimes updated by changes that users of the the local OSM server make in edit mode (by iD editor). But Overpass-API uses its own database, so I need to periodically sync these two databases.

I have already populated the Overpass-API database with openstreetmap database' data:

cd /home/rion/openstreetmap

# Dump all data from the openstreetmap database.
osmosis --read-apidb host="localhost" database="openstreetmap" \
user="rion" password="password" validateSchemaVersion="no" \
--write-xml file="main-database-dump.xml"

# Compress this dump.
bzip2 -c main-database-dump.xml > main-database-dump.xml.bz2

# Specify paths to Overpass-API's executable, database, replicate files.
EXEC_DIR="/home/rion/openstreetmap/Overpass-API/build"
DB_DIR="/home/rion/openstreetmap/Overpass-API/build/db"
REPLICATE_DIR="/home/rion/openstreetmap/Overpass-API/build/replicate"

# Upload the dump to Overpass-API database.
./Overpass-API/src/bin/init_osm3s.sh main-database-dump.xml.bz2 $DB_DIR $EXEC_DIR

So, after some changes made in the openstreetmap database, I need to upload them to Overpass-API database, too. I use this command to pull out the openstreetmap database' changes to xml file:

# Specify dates between which I will read openstreetmap database' changes.
NOW=$(date -u +"%Y-%m-%d_%H:%M:%S")
HOUR_AGO=$(date -u -d '1 hour ago' +"%Y-%m-%d_%H:%M:%S")

# Read openstreetmap database' changes between $HOUR_AGO and $NOW dates.
osmosis --read-apidb-change host="localhost" database="openstreetmap" \
user="rion" password="password" intervalBegin=$HOUR_AGO \
intervalEnd=$NOW validateSchemaVersion="no" \
--write-xml-change file="main-database-changes.xml"

# Compress the changes file.
bzip2 -c main-database-changes.xml > main-database-changes.xml.bz2

So now I have the compressed xml changes file and need to upload it to the Overpass-API database. So the question is: how to do it?

I know that there is a ./Overpass-API/build/bin/update_database util, but how to use it? I've tried to use it this way:

./Overpass-API/build/bin/update_database --db-dir=$DB_DIR \
main-database-changes.xml.bz2

But it gives me an error:

Unkown argument: main-database-changes.xml.bz2

Usage: update_database [--db-dir=DIR] [--version=VER]
[--meta|--keep-attic] [--flush_size=FLUSH_SIZE]
[--compression-method=(no|gz)] [--map-compression-method=(no|gz)]

How do I use the update_database util with a local changes file or are there some other utils to do it?


texts: After I give kring command, this error comes up: data item in gstat object and newdata have different coordinate reference systems

texts:
> library(gstat)
> library(readxl)
> library(sp)
> library(rgdal)

#2013(FATALITIES)

data <- read.csv("C:/Users/sbs/Desktop/Study/RA/Accident analysis/New folder/Iowa total crashes 2013_cleaned.csv", header = TRUE, sep =",")
plot(data)
attach(data)

> class(data)

[1] "data.frame"

> coordinates(data) <- ~XCOORD+YCOORD

> class(data)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

> variogram1<- variogram(FATALITIES~1, data)
> plot(variogram1)

> fit.variogram1<- fit.variogram(variogram1, vgm(0.002, "Exp"))

> plot(variogram1, fit.variogram1, cutoff = 300000, pch =20, col ="red")

> fit.variogram1

  model       psill    range
1   Exp 0.007896732 13527.15

Until this point everything looks good. My "data" is in CVS format, units meters. Fitted variogram model looks good.

I have made "X" my prediction surface or grid using ArcGIS having 500*500 meter as cell size of fishnet. It is a shapefile.

> x<- readOGR("H:/p/Road_Network_Info/road_fish_lable.shp")

OGR data source with driver: ESRI Shapefile
Source: "H:\p\Road_Network_Info\road_fish_lable.shp", layer: "road_fish_lable"
with 5006 features
It has 1 fields
Integer64 fields read as strings:  Id

> plot(x)

> class(x)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

Everything good until now. When I start performing kriging, this error shows up. But the plots of my "data" and Grid (X) overlap each other perfectly.

> krige(FATALITIES~1, data, x, fit.variogram1)

[1] NA
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  :
  var1 : data item in gstat object and newdata have different coordinate reference systems

I am performing analysis on fatalities in state of Iowa, USA.
Here, both X and data are in meters with same coordinates. i have checked it twice. I am using NAD 1987 Coordinates in ArcGIS to make the required shapefile.


texts:
I need to process a GPKG file in R, but it's too large (50GB) to load the entire thing into memory at once. Is there any way to incrementally load the results? Right now we are loading it with a command like this:
data <- sf::st_read(filename,
                    query = glue("SELECT * FROM layer ",
                                 "LIMIT {num_rows_per_batch} OFFSET {batch_offset};"))

in order to load only num_rows_per_batch at a time. However, when running this with a large batch_offset, it takes a long time (I'm guessing because behind the scenes SQLite is iterating through all the previous results. Rather than loading each batch with an entirely separate query, is there any way to use a single query, but load the results incrementally, so that we could process them one batch at a time? I'm thinking something like the way RSQLite::dbFetch() allows loading one batch at a time (see, ex: the "Batched queries" section in the RSQLite vignettes).


texts: How to convert GTFS stop data to a spreadsheet?

texts:
I'm trying to make a timetable in a spreadsheet by converting the agency's provided GTFS data.

Unfortunately, I have been unable to find any tools for doing this. It seems as though many people use GTFS for maps, not scheduling.

Is there any tool that can make it a spreadsheet, or even HTML?

I looked at BlinkTag's tool but I can't use node-js so it doesn't meet my requirements.


texts:
After classification of an image in GEE, I got 2 lists of user's and producer's accuracies (1*7 and 7*1 dimensions). How can I calculate an average value for each of them now?

print('Validation consumersAccuracy: ', testAccuracy.consumersAccuracy());
print('Validation producersAccuracy: ', testAccuracy.producersAccuracy());

Applying the list reducer does not work.

var meanConsumer=ee.List(testAccuracy.consumersAccuracy()).reduce(ee.Reducer.mean())

Here is the whole code (it's not actually mine, I am new at GEE)
https://code.earthengine.google.com/b9ebf87363121a6e51a980dfea187166


texts: Where can I find reflectance / VI data for crops for different dates across their growing season?

texts:
I have been trying to find literature values for reflectance of soybean for different dates along the growing season (I am specifically interested in 30, 60 and ~120 days after planting - so the V3/V4 stage, R3/R4 stage and R8 stage) but all the spectral response graphs (there's only graphs - I haven't been able to find a table or a dataset, except on one occasion) that show up in my literature search have "Average Spectral Response" and no date indicated. There is a 1984 dataset on the NASA website which has measurement dates mentioned - but no mention of planting date, and the reflectance data doesn't mention the measurement dates either.
Is there a paper or some other resource where I could find such data?


texts:
I'm trying to use QTiles to generate tiles from a shapefile layer in QGIS 2.18.1.  The tiles are for a web map, so my layer CRS is EPSG 3857.  But QGIS reports the layer extents in degrees instead of meters: xMin,yMin -111.046,32.1004 : xMax,yMax -110.777,32.2978.

In the JSON metadata produced by QTiles, the extents are extremely wrong:

{
   "maxZoom":18,
   "minZoom":0,
   "bounds":"-0.000997541174691,0.000288363206654,-0.000995126527844,0.000290135777677",
   "name":"derp",
   "format":"png"
}

It appears that QTiles expects the coordinates for an EPSG 3857 layer to be in meters (or EPSG 3857 pseudo-meters) and scales them down by a factor of approximately 111,320 to get degrees.  QTiles is rendering a tiny area around 0, 0 and producing only blank tiles.

If I change the layer CRS to EPSG 4326, QTiles produces tiles as expected.  But then isn't the projection slightly incorrect for a web map?


texts: How to programmatically import a Shapefile into Oracle DB using Java?

texts:
What's the best approach to import a user uploaded shapefile into an Oracle DB? The dbf, shx and shp files of the Shapefile come in a ZIP Archive.

Currently, we unzip programatically and then parse the dbf, shx and shp files manually using custom Java and the GeoTools library. From this we build SQL Strings to create the table and insert the values. This works, but seems cumbersome given the complexity of shapefiles.

I thought that there must be an easier way of just using some kind of Java library to which we pass the shapefile. The library should just take care of all the shapefile import stuff.

However, I do not find any GeoTools method for simply importing a Shapefile into a Oracle DB. I only found this Oracle Java Program, but no tutorials on how to use this as part of a java based web portal.


texts: Reducing filesize for landslide susceptibility map

texts:
I have a large .tif file that I would like to reduce to under 15MB (currently 650MB). The file shows 3 zones of landslide susceptibility (see pic below). I am processing it so that I can upload it to Mapbox.org, which requires smaller file sizes.

So far I have tried Boundary Clean and Majority Filter.
I then converted it to a shapefile and would like to run simplify, but there's not enough memory to do so.

I am now thinking of splitting it into sections and processing in piecewise manner. I am wondering if there's any better way to go about this. Some simplification (30-60 meters of offset) is ok, but I'd like to retain general shapes and locations of information.


texts: QGIS find and then replace value in multiple columns

texts:
This is based on a previous question/answer that I had. Replacing values in multiple columns using QGIS Field Calculator?

I need it to update 3 fields now and I get an name 'e' is not defined error. What have I done wrong?

layer = qgis.utils.iface.activeLayer()

# Set field names you want to update
field_1 = "New_Name"
field_2 = "New_File_P"
field_3 = "Process"
idx_1 = layer.fieldNameIndex( field_1 )
idx_2 = layer.fieldNameIndex( field_2 )
idx_3 = layer.fieldNameIndex( field_3 )

# Set expression to find all features which fall into expression
exp = QgsExpression( """ "Existing_F" LIKE '%2015 Aerial Ortho\\_2015\\ECW_Tiles%'  """ )

# Select all features which fall into expression
ids = [i.id() for i in layer.getFeatures(QgsFeatureRequest(exp))]
layer.setSelectedFeatures(ids)

# Set expressions to fill in values for selected features
formula_1 = """ "File_Name" """
formula_2 = """ '\\sipv-gis01\\GIS_Library\\Raster\\Aerial\\2015_25cm\\Images\\ECW' """
formula_3 = """  'y'  """
e_1 = QgsExpression(formula_1)
e_2 = QgsExpression(formula_2)
e_3 = QgsExpression(formula_3)
e.prepare(layer.pendingFields())

with edit(layer):
    # For each selected feature
    for f in layer.selectedFeatures():
        f[idx_1] = e_1.evaluate(f)
        f[idx_2] = e_2.evaluate(f)
        f[idx_3] = e_3.evaluate(f)

        layer.updateFeature(f)

Error

e.prepare(layer.pendingFields())
Traceback (most recent call last):
  File "<input>", line 1, in <module>
NameError: name 'e' is not defined


texts: Setting up data from google geocoding correctly for ST_DWithin search?

texts:
I am mostly relying on What Spatial Reference System do I store Google Map's Lat/Lng in

My web application takes user addresses, queries the Google geocoding API for coordinates, and then stores those coordinates in a postgresql postgis-enabled database. To do this I run the following insert command:

$query = "INSERT INTO userloc (email, latitude, longitude) VALUES($1, $2, $3)";
$rs = pg_query_params($con, $query, array($emailN, $latitude, $longitude));

followed by a command within postgresql:

UPDATE userloc SET the_geom = ST_SetSRID(ST_MakePoint(longitude,latitude),4326);

I am only storing data for a few hours (far less than 30 days), so my understanding is I am not violating google TOS.

What I did is based on the post I referenced saying that google geocoding already comes in 4326 format.

Once I have a collection of these points and I want to find all points within a certain range of a particular longitude, latitude pair, what query do I run?

I have tried the following:

$query = "SELECT * FROM userloc WHERE ST_DWithin(the_geom, 'POINT(".$longitude." ".$latitude.")', 100.0)";
$rs = pg_query($con, $query);

When I run this, I get the following error:

pg_query(): Query failed: ERROR: Operation on two GEOMETRIES with different SRIDs in /var/www/html/signup.php on line

If I try running this instead, using ST_SetSRID:

$query = "SELECT id FROM userloc WHERE ST_DWithin(the_geom, ST_SetSRID(ST_MakePoint(".$longitude.", ".$latitude."),4326), 100.0)";
$rs = pg_query($con, $query);

I don't get an error, but the search doesn't make sense. Or it might make sense if it I knew more about geometric units because the distance parameter (100 above) clearly doesn't correspond to normal distance. If I make it very small it does cut to the nearest (a cluster in NY), but if I make it just slightly larger it starts including points in Chicago.

I thought this probably had something to do with geography vs. geometry. I am interested in geographic searches, so then I thought I should convert to geography like so:

$query = "SELECT id FROM userloc WHERE ST_DWithin(the_geom::geography, ST_SetSRID(ST_MakePoint(".$longitude.", ".$latitude."),4326)::geography, 100000.0)";
$rs = pg_query($con, $query);

This really seems to do the trick (can distinguish nearness within cities), but the units are still baffling. I thought the default units were meters, but I had to enter the # above, 1 with 5 0's to get all the points in Chicago to be included in a Chicago-oriented search, so I don't think it is meters.

So two questions.

Can someone confirm that what I have done makes sense?

I am sorry I am so new to gis in general, so I am struggling to get a hold of all the different systems of measurements, data storage.

Second, can someone tell me what the distance is in my last command or how I can find out?

I need to be able to set a specific distance.


texts:
I am currently creating two similar maps with the same basic layers, including same basemap.  I am exporting the mxd as a pdf with the exact same parameters, yet one comes out pixelated, and one is smooth (see attachments).  I have tried different options for export from Compression to Picture Symbol (i.e. switch default from "rasterize layers with bitmap markers/fills" to "vectorize..." yet no improvement on pixelated figure.  Data is the same for figure samples.
Pixelated map is 3Mb and the smooth map is 1Mb.
Has anybody experienced a similar issue?
ArcMap 10.6

Smooth figure

Pixelated figure


texts: Geoserver SLD using PropertyIsEqualTo And PropertyIsLike

texts:
I have two SLD rules one using  PropertyIsEqualTo and another using PropertyIsLike

Using  PropertyIsEqualTo

<Rule>
          <Name>Vendor FirstVendor</Name>
          <Title>FristVendor</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>Vendor</ogc:PropertyName>
            <ogc:Literal>firstVendor</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#0033CC</CssParameter>
                </Fill>
              </Mark>
              <Size>8</Size>
            </Graphic>
          </PointSymbolizer>
     </Rule>

Using PropertyIsLike

<ogc:Filter>
            <ogc:PropertyIsLike  wildCard="*" singleChar="." escape="!">
            <ogc:PropertyName>Vendor</ogc:PropertyName>
            <ogc:Literal>firstVendor*</ogc:Literal>
            </ogc:PropertyIsLike>
</ogc:Filter>

The SLD using Like works fine where as using equal does not work at all.I have seen the string in database well there is not leading or trailing spaces.

Geoserver is taking space at the end of literal. That was the issue. I added space to PropertyName eg Vendor_


texts:
I have a countries' polygon with an attribute table. I also have a raster.
I want to multiply the values of the polygon (of each country) with the raster.

Polygon

Countries   Cost

C1            X

C2            Y

C3            Z

I also have a raster which includes these countries. Raster value is R.

I want to create another raster which will be the result of the multiplication of the polygon values with the raster.

My idea was to convert the polygon to raster and use the value I want to multiply with the raster. However, the raster comes with 0 to 255 values and not with the ones I have on the attribute table.

Working on ArcMAP 10.2, the rasterized polygon comes with 0-255 values


texts: Is there a way to add Common Fields to multiple Feature Class's in ArcCatalog?

texts: I've got a File Geodatabase that has many different feature class's.  Now versus adding common fields one feature class at a time to the attribute table is there a toolbox or way to add a field to multiple feature classes at once?


texts:
I have exported a polygon shapefile in QGIS 3.16 as a DXF by selecting 'Feature symbology' as the Symbology export options.
When I open it in Global Mapper and use 'Feature Info Tool' to select a polygon, I see the image below:

By visualizing the Feature Information, I see the attributes as below:

As shown in this image, under 'Attributes', there are two elements present: DXF-HATCH-START-LINE and ELEVATION. This file will be applied further for procedures that require it to not having any attributes and the Feature Information Attribute column has to be empty.
I am not familiar with CAD at all -- is it possible to this?


texts: Mapbox GL JS ImageSource wrong coordinates

texts:
I'm struggling to position an image correctly onto my map. I'm using one of the following URLs for my (live) image:

https://maps.dwd.de/geoserver/dwd/wms?service=WMS&version=1.3.0&request=GetMap&layers=dwd:RX-Produkt&styles=&width=1024&height=1024&srs=EPSG:3857&bbox=230596.40272276622,5934337.07299266,1750026.7710981474,7343518.710858418&format=image/png&transparent=true
https://maps.dwd.de/geoserver/dwd/wms?service=WMS&version=1.3.0&request=GetMap&layers=dwd:RX-Produkt&styles=&width=1024&height=1024&srs=CRS:84&bbox=2.0714827302884133,47.14423415016973,15.72075796095801,55.07980923136505&format=image/png&transparent=true
https://maps.dwd.de/geoserver/dwd/wms?service=WMS&version=1.3.0&request=GetMap&layers=dwd:RX-Produkt&styles=&srs=EPSG:1000001&bbox=-523.462,-4658.645,376.538,-3758.645&width=1024&height=1024&format=image/png&transparent=true

These are the associated WMS capabilities: https://maps.dwd.de/geoserver/ows?service=WMS&version=1.3.0&request=GetCapabilities

I tried adding an image source according to this example. Using these coordinates:

[2.0714827302884133, 55.07980923136505],
[15.72075796095801, 55.07980923136505],
[15.72075796095801, 47.14423415016973],
[2.0714827302884133, 47.14423415016973]

The problem is that the image is not aligned 100 % correctly, but slightly shifted/stretched. I read something about Mapbox using the Mercator projection and not being able to handle others? That's why there are three different URLs above, but all of them fail. Next I tried to convert my coordinates like this:

var westLongitude = 2.0714827302884133;
var eastLongitude = 15.72075796095801;
var southLatitude = 47.14423415016973;
var northLatitude = 55.07980923136505;

var topLeft = { lng: westLongitude, lat: northLatitude};
var mTopLeft = mapboxgl.MercatorCoordinate.fromLngLat(topLeft, 0).toLngLat();

var topRight = { lng: eastLongitude, lat: northLatitude};
var mTopRight = mapboxgl.MercatorCoordinate.fromLngLat(topRight, 0).toLngLat();

var bottomRight = { lng: eastLongitude, lat: southLatitude};
var mBottomRight = mapboxgl.MercatorCoordinate.fromLngLat(bottomRight, 0).toLngLat();

var bottomLeft = { lng: westLongitude, lat: southLatitude};
var mBottomLeft = mapboxgl.MercatorCoordinate.fromLngLat(bottomLeft, 0).toLngLat();

var mercatorCoordinates = [[mTopLeft.lng, mTopLeft.lat], [mTopRight.lng, mTopRight.lat], [mBottomRight.lng, mBottomRight.lat], [mBottomLeft.lng, mBottomLeft.lat]];

However this also didn't align the image correctly. I think, I'm using wrong coordinates. Any idea how to do it right?

Here's code to play around with (please compare the location of the radar blind spot (grey) with the image below to see the offset): js fiddle

Here's where the radar blind spot (pink) has to be:


texts: Using the Processing Framework in QGIS 2.0

texts:
I am writing a plugin for QGIS 2.0, and I want to use the processing framework to call some functions.

When I open the pyhton console inside QGIS, I am able to call it like this:

    >> import processing

And then, for instance:

    >>> processing.alglist()

However, if I try to import it in my Python Plugin, it fails to find the processing include:

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

I googled the problem and found this directives:

    from processing.core.Processing import Processing
    from processing.core.ProcessingConfig import Setting, ProcessingConfig

But unfortunately it also does not work for me :-/

Does anybody know how to include the qgis processing framework in a plugin, and how to call the algorithm (lets say for instance "processing.runalg(...)

I really appreciate your answers, as I found many explanations on how-to do it in the Python console, but not outside QGIS :-/

Thanks again for your time and have a nice day,
                                         Joana


texts: How to create filled contour lines in QGIS?

texts:
How can I create contour lines that are filled with a graduated symbol usinq QGIS?

In ArcGIS for Desktop I can do this. What about in QGIS?


texts: How to ad .ini files to QuickMapServices

texts: I'm using QGIS Las Palmas with the QuickMapServices plugin and want to add additional data sources from the contributor pack so they are available from the menu. I downloaded the pack, but how do I add the ini files to QuickMapServices?


texts: QGIS rounding values from attributes on python console

texts:
I'm trying to extract attribute data from a layer in QGIS (version 2.18.2 for mac) to use it internally on my plugin (I need that due to complex matrix operations with numpy and scipy), but QGIS is always returning rounded values with different methods, e.g.:

Original data when viewing it from open attribute table option:

[3.46000000, 0.48000000, 12.86999999, 13.30000000, ...]

Using getValues() or getDoubleValues() functions:

layer = iface.activeLayer()
print layer.getValues('G2')[0]  #G2 is the field name
>>> [3.46, 0.48, 12.87, 13.3, ...]

Using getFeatures() function:

field = [f['G2'] for f in layer.getFeatures()]
print field
>>> [3.46, 0.48, 12.87, 13.3, ...]

And using dataProvider() function:

dpro = layer.dataProvider()
feat = [f['G2'] for f in dpro.getFeatures()]
print feat
>>> [3.46, 0.48, 12.87, 13.3, ...]

At the first moment I thought it could be a display limitation from console's output (I was testing there before implementing on the plugin), but running from the plugin I had the same result.

Is there a method where I can extract the data with all the decimals?


texts:
I imported a raster into QGIS but it is not displaying correctly. The workspace projection is WGS84, Zone 39N, with units of meters. The raster is in the same projection, except in U.S. feet. The raster appears to be offset by the amount expected due to a units problem. QGIS doesn't have a compatible projection in feet, and I don't see any way to specify the horizontal units. I tried to create a custom projection, using the following specification from the projection of my raster:

PROJCS["WGS_1984_UTM_zone_39N_ft",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",1640416.666666667],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",51.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0.0],UNIT["Foot_US",0.3048006096012192]]

However, QGIS complains that this projection is identical to EPSG 32639, which is in meters, not feet.

Is there a way to specify the units or otherwise fix this issue?


texts: How to automaticaly create a buffer for a line in QGIS

texts:
I am creating some lines, and I would like to add a 500m buffer around thoses lines. But i need this buffer to automaticaly change if i modify my line around.
Is there any solution to this?


texts: Copying first record from series of feature classes into separate file using ArcPy?

texts:
I have 1000 files from which I need to copy the first record and insert into a new feature class. Additionally, is there a way to do this without copying the Field headers as well?

The end result should be a single file (I'm thinking a feature class which I would then convert to Excel) with 1000 records; every record would represent the first record of my original feature classes.

Is there any way I can automate this?

I am working on ArcMap 10.5.


texts: Using xml workspace to overwrite geodatabase in ArcGIS Server 10 with ArcSDE for SQL Server?

texts: I want to transfer a geodatabase from my development server to my client's production server. We are using SQL Server, but in different versions. Hence I cannot use Backup/restore or detach/attach. The servers don't have direct connection, so that I cannot copy the database. I used Moving a geodatabase using an XML document to export the database and then Import it to the destination. The export is done successfully, but since there are existing feature classes on the target database, I cannot import the data to the destination. I have to delete the existing features or overwrite the imported xml, but I don't know how.


texts: I am trying to test the SharpMap 1.2 libraries. I cloned the build from https://github.com/SharpMap/SharpMap and ran a Build. I encountered 56 errors under VS2015. There is virtually no info on how to build on GitHub. Is there a step I am missing?


texts: Creating summary statistics table with as many rows as iterated features?

texts:
I have a feature layer of 25 buffers, which are the light blue circles you see. Let’s call each buffer a community. Many buffers overlap. Each buffer has a unique ID.

I also have a feature layer of points (centroids). Some of the centroids fall within certain buffers, some fall outside. Every centroid has columns of information associated with it, two of which are important to me: Origins and Destinations.

WHAT I WOULD LIKE TO DO:

(1) I want to iterate through each feature in the buffer layer (loop through the buffers one-by-one) and ask: which centroids are completely within this buffer? (2) Then, I want to ask, of the centroids that fall within the buffer, what is the sum of the origins and the sum of the destinations? (3) Last, I want to combine the results of (1) and (2) in a table with three columns of information: The buffer ID, the sum of origins in that buffer, and the sum of destinations in that buffer.

MY SETUP and MY PROBLEM

SETUP

I am using the Feature Selection Iterator in ModelBuilder to loop through the buffer file. I have chosen to group by the BufferID field.
I use the Select Layer by Location tool to select the centroids that are ‘completely within’ each iterated feature in the buffer file. (The buffer is the selecting feature, the centroids are the input feature layer)
The output of this geoprocess, labeled Centroids(2) in the image below, is, I believe, the collection of centroids ‘completely within’ a particular buffer.

I use the Summary Statistics tool to calculate the sum of the Origins and Destinations associated with Centroids(2). In the Summary Statistics tool interface, I select Origins and Destinations as statistics fields and I calculate their sum. The output of this tool is a table. I select the location of this output table.

PROBLEM

This is where I am stuck. If I run the model above, I get this:

This is interesting, because I see the tools turning to red as many times as there are buffers – so I do believe that model is iterating through each feature in the buffer layer, selecting centroids by location, and creating summary statistics, when run.
Clearly, though, I have not set up the model properly to get the desired output because the resulting table has only one row. There should be as many rows as there are buffers, which in this case is 25.  Moreover, I would like the OBJECTID* in the above image to be the BufferID from the buffer layer. When I grouped the Feature Selection Iterator by BufferID I thought this might carry over to the Summary Statistics.
Although the model is iterating 25 times (as evidenced by the tools turning red 25 times) the output it generates is for one iteration only. What do I have to change in order to get (A) a table with 25 rows, and (B) an ID that corresponds to the BufferID?


texts: I am not able to display Ground Overlays in a KML file with OpenLayers. Does OpenLayers supports KML Ground Overlays? If not, how would I proceed in writing my own code or plugin for OpenLayers so that I can display Ground Overlays?


texts: Cutting TIFF in QGIS 3.0.2 produces an empty Raster

texts:
I am currently trying to cut a georeferenced TIFF image by using a Polygon as cutting boundary.

Both are saved in EPSG 25832.

When running the algorithm from the Raster - > Extract menu it runs and creates an empty Raster file, throws an error:

Error 1: cutline Transformation failed.

Also, the output Raster is set to EPSG 3044 instead of 25832.

When I open this empty file in the Foto viewer on Windows, it's just white. When I try to colour the pixels with 0 value in QGIS nothing is there. When I zoom to extend focus is on the right Spot where the Raster should be (but I have to set the right EPSG before - obviously).

How can I manage that?


texts: Geopandas spatial join result and dissolve by geometry

texts:
For my process I make a spatial join between a polygon grid and points :

join = gpd.sjoin(rotate_grid, gdf, op='contains')

It's look like this.

the result is a grid where each square of my grid is multiplied by as many points that are contained.

My goal is to have only one square and mean of a particular column, NdviM.
First I just tried to dissolve square without mean calcul with

dissol = join.dissolve(by='geometry')

But it returns me KeyError: 'geometry'

Somebody knows how to do this ?


texts: Joining points in multiple polygon features? Union with points and polygons?

texts:
I have a shapefile full of points (points.shp) and another shapefile with polygons (polygons.shp). Polygons.shp has multiple features in it, many of them overlapping with each other.

I would like the run an analysis [like a union] where every combination of points.shp and polygon.shp is returned.

Currently I am separating the features in polygons.shp into different files where there are no overlaps, since QGIS's Join Attributes By Location tool will only return a 1:1 match.

Is there a faster or better way to be doing this in QGIS?


texts: Opening layers with multiple geometry types in QGIS

texts:
I recently had some difficulties opening some MapInfo TAB vector layers in QGIS because each MapInfo layer contained two different geometry types: point and line.  The MapInfo layers had originally been converted from AutoCAD DWGs.  When opening in QGIS, only the point features appeared and not the linear features.

Is there a way for QGIS to open both geometries of a mixed-geometry layer or to convert or save them into two separate-geometry layers?

Failing either of those, is there a way to specify which geometry type should be opened (in my case, I was interested in the lines and not the points)?


texts: Pan to features drawn from $getJson request in Leaflet

texts:
I am quite new to both javascript and php.
However, I am making a website with a leaflet map letting users find pubs in a certain area.
I created a funciton called draw(), which is called when pressing the submit button after entering a text string in the search box. Here is the function:

if (map.hasLayer(postGroup)) {
    postGroup.clearLayers();
}
if (map.hasLayer(pubGroup)) {
    pubGroup.clearLayers();
}

 sokverdi=document.getElementById("search").value

$.getJSON('http://localhost/user/postcoderetr.php?search=' + sokverdi,
function(data) {
    L.geoJson(data, {
    style: function(feature) {
        str=document.getElementById("search").value;
        if (feature.properties.label==str) {
            return {color: "#ff0000", weight:1, opacity: 0.7};
            } else {
            return {color: "#33AFFF", weight:1, opacity: 0.4};
            }
        }}).addTo(postGroup);
});
  $.getJSON('http://localhost/user/pubinpostcode.php?search=' + sokverdi,
function(data) {
    L.geoJson(data, {
    pointToLayer: function (feature, latlng) {
      return L.circleMarker(latlng, markerOptions);
    }}).addTo(pubGroup);
})
;

$.getJSON('http://localhost/user/pubnamesearch.php?search=' + sokverdi,
function(data) {
    L.geoJson(data, {
    pointToLayer: function (feature, latlng) {
        return L.circleMarker(latlng, markerOptions)}
    }).addTo(pubGroup);
})};

The question is: How do I either pan to or zoom to the layer that has been retrieved from the databse via the PHP script in the $.getJSON method?
I would like Leaflet to pan to the features drawn in the map, no matter which PHP request they are coming from. Should the panning be done inside this function, or outside?


texts: installing older version of QGIS on Ubuntu

texts:
I recently installed the latest version of QGIS (2.0) on my UBUNTU 12.04. But am not very satisfied with it. I would like to uninstall it and reinstall QGIS 1.8 v on Ubuntu. Am unable to do so. It would be great if anybody could help me with this.

thanks


texts: Zones missing after running zonal statistics in ArcGIS Spatial Analyst?

texts:
I have a problem using zonal statistics in ArcGIS.

My raster contains values/data with 30m cells.

I had several vector layers, containing administrative subdivisions of Belgium. I have converted the vector layers to raster data with the Polygon to Raster tool. I used the cell center assignment and used the same cell-size than that of my value raster data.

Each time I use zonal statistics, 2 zones are missing. I don't understand why ?

It's not a size problem as these are not necessary the smallest zones.

I tried the same procedure using another raster data sets for defining the zones, and the same occurs at a different location.


texts:
I am using GeoPandas's sjoin function to join 2 dataframes: dataframeA has latitude and longitude information whereas dataframeB has polygon information. Number of rows in dataframeA may vary (~70M) but are the same for dataframeB (825k). I want to perform point in polygon operation and update dataframeA with information from dataframeB. Here is my code which works (rtree and libspatialindex has been installed):
dataframeB = gpd.GeoDataFrame(dataFromReadCSV,crs="EPSG:4326",geometry=geometry)
dataframeA = gpd.GeoDataFrame(dataframeA,crs="EPSG:4326",geometry=gpd.points_from_xy(dataframeA.longitude, dataframeA.latitude))
dataframeA = gpd.sjoin(dataframeA, dataframeB, op='within', how='left')

Since the memory requirement for this task is very high, I chunk dataFrameA before sjoin and append the results from disk. This process has been working fine.
Environment: Python 3.6; Dask - for high performance cluster
Problem: For chunked dataframeA (~7-8M rows), it takes about 2-3 hrs. I know point in polygon is computationally expensive.
Is there a way to speed this up?


texts: QGIS Desktop 3.10.9 with GRASS 7.8.3 doesn't work in Windows 10

texts:
I just installed QGIS LTR (qgis-ltr-full) from OSGeo4W in Windows 10 (version 2004), but when running QGIS with GRASS I get the following error:

qgis-ltr-bin-g7.exe - Entry Point Not Found

The procedure entry point
?glyphCache@QFontEngine@@QEBAPEAVQFontEngineGlyphCache@@PEBXW4GlyphFormat@1@AEBVQTransform@@@Z
could not be located in the dynamic link library C:\OSGEO\apps\qtS\bin\Qt5OpenGL.dll.

And after "Aceptar":

Error loading QGIS

Oops, looks like an error loading QGIS
Details:
Could not load qgis_app.dll
Windows Error: No seencontró el proceso especificado.
Help:
Check C:\OSGEO\bin\qgis-ltr-bin-g7.env for correct environment paths

What's the problem?


texts:
I want to clip raster (SRID 32649) with polygon(SRID 4326) into SRID 32648

First I transform raster using query

create view ras68 as
select st_transform(rast,32648) as rast
from ras69

Then cutting raster 32648 with polygon

select st_union(st_clip(rast,1,st_transform(geom,32648),true) )
from ras68, shp
where st_intersects(rast,st_transform(geom,32648))

The result has error:

ERROR:  rt_raster_from_two_rasters: The two rasters provided do not have the same alignment
********** Error **********

ERROR: rt_raster_from_two_rasters: The two rasters provided do not have the same alignment
SQL state: XX000

What's my mistake and how to fix that ? (If run st_intersects() only, the result are fine, they overlap)


texts:
I am getting the polygon id by doing

identifyParams.geometry = evt.mapPoint;
identifyParams.mapExtent = map.extent;
identifyTask.execute(identifyParams, function(result) {
     var i;
     var polyId= new Array();
     for(i=0; i<result.length;i++){
       polyId[i]=result[i].feature.attributes["UNIQPOLYID"];
     }

I want to zoom in to the polygon when the value of i is 1. How will I zoom in to that particular polygon?


texts: How to split string (the last digits) in QGIS?

texts:
I´m working with QGIS 2.12. Is it posible to split string "from the right side" (it means I needed the last text from every field)?

In one field I have description like:

new street 25

old street 2

lower 26A

lower new street 125

Jozefs and Elisabeths new gardens 147A

Every line has a different count of texts. I needed a column with numbers:

25

2

26A

125

147A

Could be possible in QGIS?


texts: Selecting contiguous area with raster in ArcMap?

texts:
In ArcMap 10.4, I have created a slope raster from a DEM and am trying to select only contiguous areas of slope less than 10 degrees. I have used the raster calculator to select that subset of slopes (shown here in pink), but I only want to include areas that are contiguous (i.e.   no gaps between pixels, as seen in the one linear raster area below, which is a road).

I would also like to rule out areas of the slope raster that  have significant pockets of slopes that are steeper than 10 degrees (as seen in the blobs on both sides of the road).

Is there a focal or zonal function that can do this?


texts: Adding attributes to OutTable of Check Geometry tool using ArcPy?

texts:
I am writing a script that checks the geometry of all point, lines and polygons in a gdb. I used a script provided by Esri and it works fine, however it only returns the name of the featureclass, object Id and problem.

I would also like it to return a specific attribute of each feature with a geometry problem; Namely "PUIcode"

I used the following script:

# The geodatabase in which the feature classes will be checked
source_gdb = r"C:\Users\deboerr7373\Documents\assets.gdb"
outTable = "checkGeometryResult"

# A variable that will hold the list of all the feature classes
# inside the geodatabase
def listFcsInGDB(gdb):
    ''' list all Feature Classes in a geodatabase, including inside Feature Datasets '''
    arcpy.env.workspace = source_gdb
    fcs = []
    for fds in arcpy.ListDatasets('','feature') + ['']:
        for fc in arcpy.ListFeatureClasses('', 'point' or 'line' or 'polygon' ,fds):
            fcs.append(os.path.join(fds, fc))
    return fcs

fcs = listFcsInGDB(source_gdb)

print("Running the check geometry tool on {} feature classes".format(len(fcs)))
arcpy.CheckGeometry_management(fcs, outTable)

print("{} geometry problems found, see {} for details.".format(arcpy.GetCount_management(outTable)[0],
                                                               outTable))

And I got the following results "CheckGeometryResult":

However What I would like to have:

The data in my bdg are point/line/polygon feature classes and have the following attributes:

Has anyone any idea how to solve this problem?
I considered joining tables afterwards, however the CheckGeometryResult table contains objects id's that come from different features classes. So I would than have to join on the combination of featureclass name and object Id. But I could only find the option of joining on only object Id.


texts: pgAdmin4 combines geometry SRID and type

texts:
Using PostGIS, I have a table with a geometry of LineStringZ with SRID 28356. In pgAdminIII (v1.22.2) its definition in the SQL pane shows as:

geom geometry(LineStringZ,28356)

In pgAdmin4 (v4.1) the exact same table's SQL definitions shows as:

geom geometry(7259142)

Another table with the same SRID but a PointZ type shows in pgAdminIII as:

geom geometry(PointZ,28356)

And in pgAdmin4 as:

geom geometry(7259138)

Assuming this is a feature, not a bug, it looks as though pgAdmin4 is combining the type and SRID into a single "code", e.g. (LineStringZ,28356) --> 7259142, but I can't find information on how this is performed. I have 2 questions:

How is this new "code" calculated?
Is there a way to show the type and SRID in pgAdmin4, just like it was in pgAdminIII?

Edited to add my environment:

PostgreSQL 9.6.1 on x86_64-pc-linux-gnu, compiled by gcc (Debian 4.9.2-10) 4.9.2, 64-bit
POSTGIS="2.3.1 r15264" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER
pgAdmin4 v1.1 using python 2.7.12 and Flask 0.11.1
pgAdminIII v1.22.2


texts: Set "reduce resolution" method during export in Earth Engine?

texts:
I've created an EE image object with a bunch of bands that lacks any projection information (per the GEE documentation: In particular, note that images produced by reducing an ImageCollection have no projection. This means that you should explicitly set the scale on any computations involving computed images output by an ImageCollection reduction.). I want to aggregate (that is, use a mean reducer in GEE lingo) my data to a standard, moderate resolution--say, 1km. But since .reduceResolution only accepts images that already have a projection, I find myself in the awkward position of having to force reproject to a smallish resolution, then use .reduceResolution(mean), then force reproject again to my desired size...all before having to set the scale again when I go to export the image.

Surely I am making this more complicated than it needs to be. Is there a way to set the export parameters such that it uses .reduceResolution(mean) when the data is exported?

Here's some simplified code to illustrate my problem:

    var AllData = LandsatComposite.addBands(Temp)
                                    .addBands(Elevation).toFloat()
                                    .reproject({crs: 'EPSG:32620', scale: 30})
                                    .reduceResolution({reducer: ee.Reducer.mean()})
                                    .reproject({crs: 'EPSG:32620', scale: 100});

Export.image.toDrive({
 image: AllData.select('B1', 'B2', 'B3', 'Temp', 'Elevation').updateMask(mask),
 description: "AllData",
 region: ROI,
 scale: 100,
 fileFormat: 'GeoTIFF',
 crs: 'EPSG:32620',
 maxPixels:1e13
 });


texts: Getting layers path in Layer properties window in QGIS 3.4 by PyQGIS?

texts:
I'm going to know how I can get layers full path in table of contents in QGIS3.4 when my QGIS is closed and it's running such as daemon?

Rememebr: my shapefiles are in some different folders.

This below is my code:

# If you are not inside a QGIS console you first need to import
# qgis and PyQt classes you will use in this script as shown below:
from qgis.core import *
import os
#This is is needed in the pyqgis console also
from qgis.core import (
    QgsVectorLayer
)
QgsApplication.setPrefixPath("/usr", True)
# Create a reference to the QgsApplication. Setting the
# second argument to False disables the GUI.
qgs = QgsApplication([], False)
qgs.initQgis()
layersName =[]

# Get the project instance
project = QgsProject.instance()
project.read('/media/sigma/Data/PyQGIS/pyqgis_project/pyqgis_hydro.qgz')

#print(project.fileName())

for Layers in QgsProject.instance().mapLayers().values():

    layersName.append(Layers.name())

path_to_shape_layers_in_toc =[]
for l in layersName:

    path_to_shape_layers_in_toc.append(os.path.join(QgsProject.instance().homePath(),l+".shp"))

vlayer = QgsVectorLayer(path_to_shape_layers_in_toc[2],"Countour Layers","ogr")

if not vlayer.isValid():
    print("Layers not loaded!")

I expect exact path shown in properties window for each layers in TOC, but with this code I get only layers path loaded in TOC, not full path.


texts: Map accuracy assessment by moving window in R

texts:
I would like to make an accuracy assessment using a confusion matrix between classified Landsat image and reference dataset. Maps have the same resolution and extent. I would like to evaluate the agreement by pixel-by-pixel. In many studies I have found that they used moving-kernel window (ideal 3x3) to "deal" with the Landsat pixels misregistration. However I couldn't find any approach to do use this moving window for confusion matrix evaluation in R, normally it is used to values interpolation.

Do you have any ideas how to implement moving window into classification accuracy assessment? Or am I misunderstanding this approach?

Thanks a lot,

Example:

library("raster")
clas <- raster(ncols=5, nrows=5)
values(clas)<-c(2,1,0,0,1,
         2,1,0,0,0,
         1,0,0,0,0,
         0,0,0,0,0,
         2,1,0,0,2)
reference <- raster(ncols=5, nrows=5)
values(reference)<-c(2,1,0,0,0,
         2,1,0,0,1,
         0,0,0,0,0,
         1,2,0,0,2,
         0,0,0,0,0)

EDIT:

After utilisation of movingFun (raster) I think this is mostly to get information in one cell taking into consideration the cells around.

so if I have a

ref<-c(0,0,0,0,1,0,0,0,0)
class<-c(0,0,0,0,1,1,1,0,1)

and I apply movingFun

mov_ref<-movingFun(ref, n=3)
mov_clas<-movingFun(class, n=3)

I will obtain

mov_ref
NA, 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333 0.0000000
0.0000000        NA

mov_clas
NA 0.0000000 0.0000000 0.3333333 0.6666667 1.0000000 0.6666667
0.6666667        NA

For accuracy assessment needs could I directly compare mov_ref and mov_class? Or do I have maybe to reclass it into 0,1 (by threshold 0.5) to compare them? If I leave the numbers as 0.3333 I dont reach a good result in classification accuracy.


texts: Rendering maps in graphic design software

texts:
Is it possible to see a map in adobe illustrator with zoom rendering? What I want to do is use google maps or OSM or any other map as a reference background for a map I am trying to create in illustrator.

I have found vector maps from openvectormaps.com, but these are still vector images-the only difference between a raster image and this is it doesn't get pixelated.. I want to render maps as they zoom, like we can do in GIS software.

I will ask this on the graphicdesign stackexchange as well.


texts: Giving specific options to user when creating tool in ModelBuilder?

texts:
I am making  a terrain analysis tool for military use (educational purposes). I am doing a raster analysis, using standardised terrain data in combination with a surface raster to determine the maneuverability in a user specified area.

My model gives a weighted raster output, which I will symbolise in a logical way if/when I get the results I want. Over time I will add more input (especially soil type) to the model, but this is the general design of it.

I am struggling to find a way to let the user differentiate between wheeled and tracked vehicles, and dry and wet terrain. My assessment is that these two aspects will affect the result the most, but so far my solution has been to make several different tools with different weighting.

My best guess would be that I have to let the user choose between these options in the weighted overlay:

Wheeled vehicles, dry soil
Wheeled vehicles, wet soil
Tracked vehicles, dry soil
Tracked vehicles, wet soil

Is this possible in any way? Or any other workaround that removes the need for several tools?


texts: QGIS QChainage doesn't create point at intervals as I input

texts:
I want to create points at 10 meters intervals along the road, but when I input 10 meters into QChainage in QGIS as follows:

I got the points at 8.33333 meters intervals. I don't know why and how to create points at intervals that I want.

here is the information of my roadnet:


texts: A way to generate random/fake sensor data

texts:
First, what I'm trying to achieve is to implement a gis-based platform in the context of a smartcity. So basically I'll need many layers, some are static like gaz,water circuits, hospitals ...., some are crowd-sourcing based (citizen's reports). The third type of layers I need is sensor-based.

Now for the last layer my problem is that I don't have actual sensors at hand. So what I'm looking for is : is there a sensor emulator, or a sensor data generator of some kind that could feed me random/fake data (temperature,CO2 pourcentage ...) ?

I'm using geoserver as my back-end and I plan to communicate with it using the WFS protocol using openlayers/leaflet as a client. And I can imagine the existance of a program that periodically sends data to geoserveur (via WFT-T maybe, or direct access to the PostGIS database). Which would simulate the actions of a sensor. I read about a protocol called SOS, I don't know if I need it in my case.

Notes :
- This is my Master's final project and I have about 2.5 months left
- My tool choice is not final (In case another backend/frontend would help me better)
- For now I'm looking for sensors in fixed locations, which will be displayed as points or markers on the map.

So to sum up, I'm seeking guidance on how to simulate a sensor and feed the results to a backend GIS server.
Any advice is welcome


texts: Where is the result of QGis Zonal Statistics?

texts:
I just ran the Zonal Statistics in QGis, the progress bar went up to 100%, disappeared, and nothing happened. Where is the result? Or does this mean it crashed and there was an error? I haven't seen any new messages in the log (View > Panels > Log Messages).
EDIT: I found it! It has been added to the polygon layer. But this is horrible!!! I don't want it to modify my polygon layer!!! I want the data separate! How can it even think of modifying my layer without asking?
Is there any way to let the Zonal Statistics put the data separately?


texts: Requesting a JSON format from GeoServer using Openlayers and GeoExt

texts:
I am requesting a vector layer from GeoServer using WFS request in JSON format. I have this URL http://myip.com:8001/geoserver/workspace/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=workspace:mylayer%20layers&maxFeatures=50&outputFormat=application/json.

I want to populate my grid/table in my client side bind to the vector layer. What happens is that I won't be able to see my layers. My question is that, am I still going to use a server proxy?

I am using Tomcat 6 in Windows 7.

Try to see my code:

var vector = new OpenLayers.Layer.Vector("My Layer");
store = new GeoExt.data.FeatureStore({
    layer: vector,
    fields: [
        {name: 'name', type: 'string'},
        {name: 'description', type: 'string'}
    ],
    proxy: new GeoExt.data.ProtocolProxy({
        protocol: new OpenLayers.Protocol.HTTP({
            url: "http://myIP:8001/geoserver/workspace/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=workspace:mylayer%20layers&maxFeatures=50&outputFormat=application/json",
            format: new OpenLayers.Format.GeoJSON()
        })
    }),
    autoLoad: true
});

var gridPanel = new Ext.grid.GridPanel({
    //title: "Feature Grid",
    region: "east",
    store: store,
    width: 'auto',
    height: 100,
    columns: [{
        header: "Name",
        width: 200,
        dataIndex: "name"
    }, {
        header: "Description",
        width: 100,
        dataIndex: "description"
    }],
    sm: new GeoExt.grid.FeatureSelectionModel()
});