texts: Problems defining custom CRS in QGIS - Orthographic projection

texts:
Working with QGIS 3.14 on Win 10, I try to define a custom CRS - I want to create an Orthographic projection centered on the Mediterranean - 35 N 20 E.
What I did: I copied the WKT-definition of North Pole Orthographic ESRI:102035 from the project properties CRS dialogue and than used this definition as starting point, pasted it in the custom CRS dialogue and changed the values from
PARAMETER["Latitude of natural origin",90,

and
PARAMETER["Longitude of natural origin",0,

to
PARAMETER["Latitude of natural origin",35,

and
PARAMETER["Longitude of natural origin",20,

So the changed WKT-definition looks like this:
PROJCRS["Mediterranean_Orthographic",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["North_Pole_Orthographic",
        METHOD["Orthographic",
            ID["EPSG",9840]],
        PARAMETER["Latitude of natural origin",35,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",20,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - north of 0°N"],
        BBOX[0,-180,90,180]],
    ID["ESRI",102035]]
Proj4
+proj=ortho +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
Ausdehnung
-180.00, 0.00, 180.00, 90.00

Clicking the validate button, it tells me the defintion is valid. When I try to click OK to save the changes, an error message pops up:
PROJ-string ist equivalent to ESRI:102035
Please try to change the CRS-definition in WKT-format

Found this here, but it did not really help: Create new custom projections in QGIS 3.12?
How can I keep the original North_Pole_Orthographic (ESRI:102035) definition and add an additional one, just centered on another point?


texts: different way to build qgis python standalone application

texts:
There seems to be 2 different ways to start a Python QGIS standalone script.
The official tutorial is using :

from qgis.core import *
# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/path/to/qgis/installation", True)
# load providers
QgsApplication.initQgis()

But I have also seen code looking like this :

from qgis.core.contextmanagers import qgisapp
with qgisapp() as app:
    #do stuff

For example here and here.

I didn't see any reference to the latter method on the official doc.

Does anybody know if there is a difference between the two ?


texts:
I have a webADF based website that I am trying to deploy. It's giving me this error. I have made sure that I have installed all the components that come from ArcGIS Server Installation media (GISServices, WebADF Runtime, WebApplications)

My server environment is: Windows Server 2008 R2, ArcGIS Server 10, IIS7

I am not sure what else it needs to find the DLL.

Apprecaite any feedback.
Thanks
Jay


texts: How to determine if two points are on the same side of a LineString?

texts:
I have two points, A and B, and a line string L. A line string is formed by a sequence of points such that consecutive points are joined by a straight line. This is the same as a LineString in Geos - I don't know if CGAL has an equivalent construct.

My question is - what is the most efficient way to determine if A and B lie on the same or different side of the polyline L? I would ideally like to use the Geos::geom library, since my data is in that form, but a CGAL solution would also work.


texts: Google Earth Engine: eq filter not taking values in list

texts:
I want to filter a raster in GEE with multiple values. I input a list of values and map over the function to apply the filter. Here is the GEE code to loop over a list of numbers and filter the raster (DEM) for values equal to the numbers in the list.

var aec = function(n){
  var ii = ee.Image(n)    // casting the number to image before filtering but didn't work
  var DEM141 = DEM.eq(ii)

  var DEM141Count = DEM141.reduceRegion({
    geometry: ROI,
    scale: 30,
    reducer: ee.Reducer.sum()
  })
  return DEM141Count
}

var elevs = ee.List.sequence(100,150,1);
var areas = elevs.map(aec)
print(areas)

I keep getting this error:
Image.lte, argument 'image2': Invalid type. Expected: Image<unknown bands>. Actual: Float.

The GEE code link is: https://code.earthengine.google.com/6e55360bc82695f930d64ee70a415dc0


texts: Losing 3d visualization of DEM in ArcScene?

texts: I am trying to see 3D of DEM using arcscene but when I change from no elevation value from surface to floating on surface dem disappear from screen does not appear even after clicking on zoom to layer?


texts: Batch export to shapefile by attribute field in QGIS

texts:
I have a project in QGIS where I am creating a map of plants over a large area. I want to take the "Genus" field (of which there are 100 hundred values) and essentially create 100 shapefiles based on each individual genus name. Eventually I'll have it set up in Geoserver and set filters from a main postGIS database, but for now I'm just playing with shapefiles for concept.

Are there any efficient ways to do this through QGIS, or are there other options that will accomplish the same purpose?


texts: Is it possible for 2 or more PostGIS extensions to coexist?

texts:
When we issue CREATE EXTENSION postgis on our PostgreSQL database, the extension created is the one of which we have the appropriate files of in our PostgreSQL server configuration.

Well, is it possible to create a situation in which we can install both postgis 2.3 and postgis 2.4 extensions for example? Maybe we can't have both extensions active at the same but perhaps we can still switch between 2 versions effortlessly, such as issuing some kind of CREATE EXTENSION postgis2.3 and CREATE EXTENSION postgis2.4 command? I am sure developers could use this to see if a PostGIS feature of a particular version is also functioning the same on another version as an example.


texts: Leaflet Event Craziness: jQuery click events don't propagate to L.DomEvents

texts:
I noticed that calling .click() on jQuery selected DOM elements in Leaflet ( such as anchor icons ) will not do anything -- the click event isn't propagated to the events created with L.DomEvents. For example, go over to bboxfinder and try to trigger the zoom out button's click event:

$('a[title="Zoom out"]').click();

Nothing will happen, the map won't zoom out. I jumped down the rabbit hole to figure out the why(s) and ended up with more questions than answers. My two favorite questions are listed below.

I'm hoping the answers for each question address two parts:

Why is this behavior happening? Be as technical as you can so people can learn!
What are the benefits of this behavior? This is a well-written, abstracted library. I assume they made these choices for reasons that are not easy for someone like me to see. Any guesses what those benefits are? For example, maybe these choices work well with mobile devices.

Questions

1. Why are the jQuery events not propagating to DOM elements? I can call click on the unwrapped jQuery DOM element ( see below ) and fire the event, but i assume jQuery would pass this event down for me:

// this works fine, notice we're grabbing the raw element
$('a[title="Zoom out"]')[0].click();

2.  When i was digging around in the Leaflet source I noticed button events were created/bound with this code workflow below ( it also happens in 3rd party packages like L.Draw ). The question is why are we using all these events on a single button and how are they working together? And really, why is it like this? I've pulled the code below from this part of the source:

    L.DomEvent
     .on(link, 'click', stop)
     .on(link, 'mousedown', stop)
     .on(link, 'dblclick', stop)
     .on(link, 'click', L.DomEvent.preventDefault)
     .on(link, 'click', fn, context)
     .on(link, 'click', this._refocusOnMap, context);


texts: Accessing elements Linestring and points from a DGN file and displaying both layers on QgsMapCanvas

texts:
I am using PyQGIS APIs to develop a standalone application which currently display's map data layers on QgsMapCanvas from .shp, .tif, .dgn files.
I am able to do it for .shp and .tif but not for  .dgn
Basically, my .dgn file contains some cities marked as points and some boundaries when i import the .dgn file in QGIS it asks me to select two elements, that are- LineStrings(boundaries) and Points(cities)
.dgn file when opened in my standalone application can only load one element i.e the boundaries(Linestrings).
How can i resolve this issue programatically using PyQGIS? Also, what type of format is .dgn raster or vector(I know it's a CAD file but it oughta' be one of raster/vector too).


texts: GeoTiff exported using GRASS will not open in other programs

texts: I am trying to export an image group as a multi-band GeoTiff using the r.out.gdal module in GRASS. For some reason I cannot get ArcGIS to recognize my GeoTiff exports. I have tried opening the file in Opticks and SAGA and neither of them can open it either. To simplify things, I generated a classification raster and tried to export this, with the same result that nothing recognizes it. I have tried playing with the module settings and adding createopt strings as recommended in the GRASS documentation for r.out.gdal, but nothing is working. Has anyone encountered this before?


texts:
I have imported a shapefile that contains a serial date number that looks like this: 41256.260799. I would like to create a new field that shows only the hour and second like this: 06:33. I've tried the available VBS and Python date functions, but nothing seems to work.

Would I make my new field of type text?


texts: GDS Workspace Path resovles as owner not user

texts:
I am having trouble retrieving the correct path to the current map layer dataSource path. I am using ArcGIS Desktop 10.2 (10.2.0.3348)

The workspace path that gets returned is correct apart from the current users name is not used but replaced by the owner of the Geodatabase Server. Which does not exist on the users machine and so returns empty dictionary on commands like arcpy.da.ListSubtypes(lyr.dataSource).

e.g.

mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, '*.Couplings'):
    print lyr.dataSource

Which returns;

C:\Users\<owner>\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\GWMAP_SQLEXPRESS.gds

What I expected would be;

C:\Users\<user>\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\GWMAP_SQLEXPRESS.gds

The user owner is the owner of the Geodatabase Server as listed in the Properties->Administration tab.

I have tried setting the default environment settings in ArcMap & ArcCatalog. Plus setting the arcpy.env.workspace to the correct path.

What am I doing wrong? Obviously I could easily use string replace to correct this but would prefer if there is a better approach.


texts:
I have US National Atlas (EPSG:2163) file having resolution of 250 m, 250 m that I want to reproject to EPSG:4326 with the same pixel size as the original map. Using

gdalwarp -overwrite -s_srs EPSG:2163 -t_srs EPSG:4326 -of Gtiff "projected.tif" "reprojected.tif",

converts the pixel size to '0.00436162,-0.00436162' which is much more area than USNatAtlas pixel area. Why is this error taking place and can I have a work around for this problem?


texts: Concatenating fields into a new field-QGIS

texts:
I'm not sure what I'm missing here. I'm trying to use the field calculator to concatenate multiple fields into one new. All fields are identified as strings.

My expression:
concatenate(  "Building"  || '-' ||   "Floor"  || '-' || "Area2" || '-' ||  "Room Numbe"  )

Expected Results:  6-01-A-02

Actual results:
6-01-A-026-01-A-016-01-A-096-01-A-106-01-A-086-01-A-066-01-A-05

Basically the calculator is concatenating all records across the above fields until reaching the character limit. I'm not sure how to only concatenate by record.


texts: Query WMS for available layers using Openlayers

texts:
Is it possible to use Openlayers to return an array of available Layers on a WMS?

I am attempting to allow users to add a layer.wms to a map.  I would like to allow the user to view available layers, select those layers that they want, and then pass that arrary to layer.wms params, add the layer and redraw.

Just wondering if there was a simple way to query layers available from a WMS.


texts: Can I use a polygon to show a WMS-layer inside while maintaining transparency outside?

texts:
I am making a map where I want a polygon filled with data from a WMS-service. What is outside the polygon from the WMS should be transparent. I am using QGIS 2.18.4.

It is a bit similar to a question on How do I clip OSM basemap with a polygon? here on GIS Stackexchange, but the solution there (inverting a polygon manually) is not sufficient. Neither would the more up to date visualization tool "Inverted polygons" be of any immediate help.

Example 1) In this example I have two WMS-layers. One in colors and one in gray. I want to use the polygon to separate the gray and color WMS-layers so that the focus area is in colors while the area around is gray. It can be accomplished in part with setting the color to gray and then playing around with transparency and/or Layer blending mode. You get the idea, but the other (gray wms-layer) does not get into play.

Example 2) This is the clearest illustration. Here I have a marine management areas layer where I want to highlight fisheries activities. The fisheries activities are from a WMS-layer with the Norwegian Coastal Administration. I do not want to show the fisheries activities in the area outside the marine management areas.

Using the first example the below image shows how I am stuck at where I get to mask out the area around the focus area. I want the blue area to mask out the color wms, and only show the gray one.

I can not find any way of connecting the masking layer to just one of the layers underneath. Is there a way to do this?


texts: Multi copy / paste point shapefile with fixed offset

texts:
I'm looking for a tool in QGIS to copy a point location of a shapefile for multiple distances.

So, for example, the point location is on X/Y coordinate, 0,0. I want to have a copy of 4 extra points on coordinate: 0,10, 0,20, 0,30,0,40 in one time.

Is that possible in QGIS 2.14 or QGIS 2.18. I'm at the moment a 2.14 user. And is there a plugin to solve this problem, I did not find a good one (https://plugins.qgis.org/plugins)?

I've looked at similar issues on gis.stackexchange but found only of offset for a single line (no multicopy).
Creating parallel lines at particular distance from lines which are part of a shapefile


texts: Finding the percentage overlap between a vector file and GPS points in QGIS

texts:
I am trying to find a way to find the percentage overlap between a vector polygon map and GPS points.
The purpose of this is to understand area usage of the animal as well as resource selection.
The points are gathered over 2 weeks and the GPS sends a signal every 2minutes when active.
My vector file is an area map that shows different area types that are of interest to my assignment. What I want is to have a table showing percentage of area used highest to lowest based on point location, together with my map.


texts: How to transform X and Y coordinate on a delimited text imported layer

texts:
I added excel sheets into my project using the deliminted tab and used the X and Y columns that are in a mercator coordinate system, I am using QGIS 3.14. These layers are not showing up in my project because another layer that I have in the project is in a different coordinate system. Showed in the picture below. The attribute table for this layer has the same X and Y coordinates as the excel sheets I brought into the project.
How do I transform the layers I brought in to match the transformation done to the layer that was supplied to me?

The layers that I imported from excel and the layer that was supplied to me have the same X and Y coordinates in the attribute table. Would saving the attribute table from the supplied and bringing it back into the project through the delimited tab work to fix the transformation issues?


texts: Convert CSV with a column that has a collection of coordinates into a spatial file

texts:
I have a CSV that has the Jsonata ( $string($eval(Location).coordinates) ) representation of the feature in a column. It is not stored in a .json file.
Can we import this into QGIS so that it has the attribution and spatial representation?
1   [Line]  [[145.44718131423,-17.2122806274458],[145.447615496814,-17.2116429871501],[145.448055714369,-17.2121269022115]] Test Line
2   [Polygon]       [[[145.443768203259,-17.2042553589467],[145.447121970356,-17.202962405607],[145.447320118546,-17.2054797649401],[145.443768203259,-17.2042553589467]]]  Test Polygon

I have tried to import using Delimited Text:WKT but get an error as this isn't valid WKT.


texts:
I'm trying to capture MouseClick events on the QgsMapCanvas. The goal is to intercept the right mouse button and display a custom context menu.
If the right mouse button is pressed with a modifier key (SHIFT+CTRL for example), the active tool action should be triggered.

I tried to install an event filter on the QgsMapCanvas provided by the iface, but this is not working.

class MyPlugin:
    def __init__(self, iface: QgisInterface):
        self.iface = iface
        self.eventFilter: MyEventFilter = None

    def initGui(self):
        self.registerEventFilter()

    def unload(self):
        self.unregisterEventFilter()

    def registerEventFilter(self):
        if not self.eventFilter:
            self.eventFilter = MyEventFilter(self.iface.mapCanvas())
        self.iface.mapCanvas().installEventFilter(self.eventFilter)

    def unregisterEventFilter(self):
        if self.eventFilter:
            self.eventFilter.deleteLater()
            self.eventFilter = None

and the EventHandler looks like this:

class MyEventFilter(QObject):

    def __init__(self, parent: Optional[QObject] = None) -> None:
        super().__init__(parent)

    def eventFilter(self, caller: QObject, event: QEvent) -> bool:
        if event.type() == QEvent.MouseButtonRelease:
            QMessageBox.information(None, 'MyEventFilter', f'event: {event}')
            return True
        return False

But when I click, nothing happens. I guess the active tool itself is already filtering the click?

I've seen that QgsMapCanvas has signals for the mouse position as well as keypresses, but unfortunately no mouse click signals.

I know that you can implement a map tool, but the idea is to have a global contextmenu handler that is always called first. Only with a modifier key the original tool action should be triggered.

Is there any mechanism to achieve this?


texts:
I have created an option in my Leaflet.js web map that exports clicked lines as GPX, but it is just exports only the last clicked line.
I need to export all clicked/highlighted/selected lines in one GPX file.

Here is my code:

$.getJSON('data/routes/Download/Nwn.geojson', function (collection) {
    Nwn = L.geoJson(collection, {
      'style': function () {
        return {
          'color': '#FFFF00',
          'weight': 3,
          'opacity': 0
        }
      }, onEachFeature: function(feature, layer) {
        layer.on('click', function() {
            document.getElementById('exportHiking').onclick = function() {
                var toJson = layer.toGeoJSON();
                var dGpx = togpx(toJson);
                document.getElementById('exportHiking').setAttribute('href', 'data:text/json;charset=utf-8,' + dGpx);
                document.getElementById('exportHiking').setAttribute('download','route.gpx');
            };
        })
        var timesClicked = 0;
        layer.on('click', function() {
          Nwn.resetStyle(layer);
          layer.setStyle({color: '#FFFF00', opacity: 1});
          timesClicked++;
          if(timesClicked>1) {
            Nwn.resetStyle(layer);
            timesClicked = 0;
          }
        });
        //   Nwn.resetStyle(layer);
      },
      distanceMarkers: { showAll: 11, offset: 1000000, cssClass: 'some-other-class', iconSize: [16, 16] }
    }).addTo(NwnGroupe);
    map.addLayer(NwnGroupe);
});

How it could be done?


texts: I'm working with QGIS and have a big DEM raster layer with height information and a vector layer with many polygons, representing 1 x 1 degree tiles. Now I want to find for every polygon the highest value of the DEM pixels in it.


texts: I need to know the value of the sum of all pixels in a raster. However, in the Statistics section of the Metadata tab in Properties, under Sum of All Cells it always indicates 0.00. I have tried with different raster formats to no avail. What am I doing wrong? How can I get this value?


texts: How to: script convertion from SRTM to 3D printable file?

texts:
Given a STRM3.tiff file as input.

How to convert this STRM.tiff file into a 3D printable file ? which I may load into most 3D printers' sofwares.

My objective is to print some terran, lunar, or martian reliefs with a basic 3D printer.

Gdal approach, QGis, or other welcome. Scriptic / terminal approach favored.


texts:
I am getting my Localisation on Leaflet map using the extension. Leaflet.Locate

When I am on my local server it works fine like shown bellow:

But when I am on a remote server I get this error.

How to solve this.


texts: Getting the output layer reference returned by Processing tool

texts:
After running Processing algorithm ending up with a layer by runalg or runandload in QGIS2, run or runAndLoadResults in QGIS 3, these methods return a result which is mostly a dictionary like {'OUTPUT': 'file_path'} etc. So, after result = processing.run*(...), I can get the output by result['OUTPUT']. But its value is mostly just a string.

If it would be a method like processing.run*(...) that returns an instance of output instead of file path, we would get the layer reference by writing just layer = processing.run*(...).

Is there any processing.run*() method that returns layer instance/reference directly, e.g. QgsVectorLayer?


texts: Why did my shapefile distance get smaller after changing the projection?

texts:
Possible Duplicate:
  Layers having the same coordinate system should overlap, but they don't. What's wrong?

I'm working with a shapefile that has a geographic coordinate system and I need to change it to a projected coordinate system to use spatial analysis tools.

Steps:

1) I defined the coordinate system using the define projection tool to State Plane (US Feet) California Zone 6.

2) Then I used the Project tool to make define the coordinate system as the state plane system using a different name. I made sure not to use the same name and followed the guidelines listed here: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00130000004v000000.htm

When I measured the county with the measurement tool, it said it was less than a mile across (it measured about 34 miles across when the coordinate system was defined using the GCS). I'm not sure what I'm doing wrong.

Any help is greatly appreciated. Thanks!


texts:
I'm developing map with Mapserver to be rendered in tiles and testing output with Ka-map. Most of the time labels split on two tiles are matching, but sometimes they're not.
I configured Ka-Map like

$tileWidth = 256;
$tileHeight =256;
$metaWidth = 4;
$metaHeight = 4;
/* $metaBuffer = Buffer size in pixels to add around metatiles to avoid
 * rendering issues along the edge of the map image
 */
$metaBuffer = 256;

but even large buffer of 256 pixels doesn't protect label from occasional splitting.
What good practices do exist to make sure labels are placed deterministically?


texts:
We have an older running ArcGIS 9.3 asp.NET application. Now we have installed 10.1 and started migrating. We had used Map Web Controls for displaying map in asp.net.Because there are some libraries deprecated in 10.1, in 10.1 environment, we got so many compilation errors. We solved one-by-one by commenting those codes. We do following things

Removed all Map Web Controls/Toolbar Controls and other controls.
We made a dummy class as "Map" with same properties and method for simulating older code and removing compilation errors.
Removed "using ESRI.ArcGIS.ADF.Web.UI.WebControls" references.

My question is ,

How much effort will required for migration?
How can we compare with the older system, because 9.3 is no more available?
We have used so many "ESRI.ArcGIS.ADF"'s client side apis also. What are the equivalent of those functions in 10.1 javascript apis?
Is any good documentation available related to 9.3 to 10.1 migration.


texts:
I am trying to create an isochrone that would show all the points that can be reached within a certain distance time, starting from a central point. I have installed pgRoutine (v. 2.2) as per instructions.

Now I'd like to install the driving_distance command as suggested here but the required files are not available anymore on Github.

Any suggestions? Or recommendations on other ways to proceed?


texts: In QGIS is there a keyboard shortcut to open/close the layers panel?

texts:
Is there a keyboard shortcut, or even something like a toolbar button, which would make it easier and quicker to open (and close) the layers panel in QGIS (2.6).

I don't see one listed under the keyboard shortcuts list.


texts:
I have two shapefiles

Blue network
yellow locations

as shown in the below figure.
I want to connect all yellow dots with the blue lines blue as shown in red color in the figure.
Is there a way to connect all the yellow circles (buildings) to blues lines(roads)?


texts:
I have to create some species distribution maps, but I have very little experience with ArcGIS unfortunately.

I have two different kinds of excel sheets:
 1. A version with many thousand X and Y coordinates
 2. A version that is identical like 1), but has an additional third column saying how many species occur at every X/Y coordinate.

I was told the first one would be the right data sheet, but don't know it that is true.

Anyway, what I need to do, is to somehow import these data into ArcMap and create a species richness map, showing where species occur and with a color-code, how many there are. That is it. This does not sound too complicated, but I can't get it to work.

Could somebody help me and tell me what to do?


texts: Merging polygons of regions within SP DF and multiple SP DFs with different levels using R?

texts:
This is actually a two-part question. I'm looking to merge some regional polygons together in R:

library(raster)
Finland_map <- getData(name = "GADM", path = tempdir(), country = "Finland", level = 2)
Aland_map <- getData(name = "GADM", path = tempdir(), country = "Åland", level = 0)

One merge I'd like to do is to include Aland as a region into the rest of Finland. Unfortunately, I can't use rbind(Finland_map, Aland_map) because the levels of the files I'm using are not the same. How do I get past that?

The second merge involves a couple of regions within Finland. Eastern Uusimaa and Uusimaa were amalgamated in 2011. I've seen similar questions and answers on the topic (such as this one), but I'm not sure on how to use raster::union() when the two polygons are within the same SP DF.


texts: WMS with additional parameter in GeoServer

texts:
In my database I have table which contains data for many users (identifying by USER_ID attribute) and I am looking for a way to filter data in my WMS service by this attribute.
For example request as below:

localhost:8080/geoserver/my_location/wms?
    service=WMS
    &version=1.1.0
    &request=GetMap
    &layers=my_location:my_location_history
    &styles=
    &bbox=54.406021,18.531216,54.497708,18.577851
    &width=768
    &height=390
    &srs=EPSG:4326
    &format=application/openlayers
    &user_id = 1

should return only data where in my tabele USER_ID = 1.

Currently, for this purpose, I am using ELEVATION dimension. If you have better and/or "more professional" solution, let me know ;)


texts: QGIS is filtering out the data I want to export

texts:
I am working on a project involving a lot of data, which I have to export from .dbf to CSV, but QGIS fails to do so and only exports less than half of the data.
It seems like QGIS is applying some internal filter that I cannot get rid of, as it shows that the data is filtered even when Show All Features is selected.

How could I solve this?
I've already tried using python with dbfread and simpledbf, but those failed to even open the file (I guess it's because they're trying to use dBase V instead of dBase III, which this file is based on).
P.s. I am using only Linux on my system, but if there are any solutions on Windows, I would be willing to give it a try.


texts:
I am trying to remove a couple elements from my webmap legend (the rest service title and the world map).

Here is what my legend code looks like:

var legend = new Legend({
     map: map,
     layer:D_ml.visibleLayers
}, "legendDiv");
legend.startup();

This is how I set up my D_ml.visible layers :

C_ml = new ArcGISDynamicMapServiceLayer("http://imcltpdarcp01.duke-energy.com/ArcGIS/rest/services/TEST_DEC_DEI/MapServer");//,{opacity: 0.6};

D_ml = new ArcGISDynamicMapServiceLayer("http://imcltpdarcp01.duke-energy.com/ArcGIS/rest/services/TEST_DEC_DEI/MapServer",{opacity: 0.75});

map.addLayers([C_ml,D_ml]);

C_ml.setVisibleLayers([2,3]);
D_ml.setVisibleLayers([0,1]);


texts: Hook for deinstallation of a QGIS plugin

texts:
I have a plugin that creates a lot of settings in the Options | Advanced dialog.

I'd like to clean up a little on plugin install.

Is there a hook I could use ?

unload() is run on any QGIS shutdown and therefore not feasible.

Thanx

/B


texts: Viewshed with RpyGeo "rpygeo.Viewshed.sa"

texts:
This function is not working:

rpygeo.Viewshed.sa(mnt, pt,out.raster="viewshed")

*Error in as.character.default(<S4 object of class "RasterLayer">) :
  no method for coercing this S4 class to a vector**

Could someone help me on this?


texts: Techniques for reducing a MapInfo TAB's file size without losing too much integrity?

texts: What are some techniques for reducing a MapInfo TAB's file size without losing too much integrity? The file needs to be reduced for use, not sharing (i.e. zipping).


texts: Find distances between points through a raster

texts:
I have a set of 79 points located overlaid on a raster that defines distance through coastal waters:  each cell equals 100 metres, with NoData for coast and islands.

I can use this raster to create a new raster with the distance of each cell to the nearest point using both the Cost Distance or Path Distance tools.

I would like, however, to identify for each point the identity of the nearest other point and the distance to it through the raster.

I am using ArcMap 10.0 with access to the Spatial Analyst tools.


texts: gdal_translate -projwin (trying to clip) creates bad file with NAN data values

texts:
I have had no luck creating rasters that are clipped from larger ones using gdal_translate.  If I just do a simple conversion (no clipping) to .tif there is no problem  However if I use -projwin and specify a smaller region I get a useless raster with naNs.

What is even more frustrating is that some raster files do work ok, especially if they are rgb.  The raster files I have to clip/convert are single band Float32 grayscale.

I am running QGIS on my Windows 10PC with QGIS v2.18.13 and GDAL/OGR 2.2.2

Here is the actual command:

gdal_translate -projwin 711541.57665 20664.3625 712593.227737 20044.0065217 -of GTiff "D:/Composite Project/testin.tif" "D:/Composite Project/cliptest.tif"

Likewise the info on this particular input file is below.  BTW - this problem is not confined to just this file but a large collection of grayscale tif, cub and img files I have:

Driver: GTiff/GeoTIFF
Files: D:/Dropbox/Composite Project/testin.tif
       D:/Dropbox/Composite Project/testin.tif.aux.xml
Size is 10032, 55134
Coordinate System is:
PROJCS["EQUIRECTANGULAR MOON",
    GEOGCS["GCS_MOON",
    DATUM["D_MOON",
        SPHEROID["MOON_localRadius",1737400,0]],
    PRIMEM["Reference_Meridian",0],
    UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",23],
    PARAMETER["standard_parallel_1",0.7],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (13140.399999999999636,27728.479999999999563)
Pixel Size = (0.520000000000000,-0.520000000000000)
Metadata:
  AREA_OR_POINT=Area
  DATA_SET_ID="LRO-L-LROC-5-RDR-V1.0"
  INSTRUMENT_ID="LROC"
  INSTRUMENT_NAME="LUNAR RECONNAISSANCE ORBITER CAMERA"
  PRODUCER_INSTITUTION_NAME="ARIZONA STATE UNIVERSITY"
  PRODUCT_CREATION_TIME=2013-05-30T01:32:18
  PRODUCT_ID="NAC_ROI_APOLLO11LOA_E004N0235"
  PRODUCT_TYPE="RDR"
  TARGET_NAME=MOON
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   13140.400,   27728.480) ( 23d26' 0.15"E,  0d54'51.94"N)
Lower Left  (   13140.400,    -941.200) ( 23d26' 0.15"E,  0d 1'51.74"S)
Upper Right (   18357.040,   27728.480) ( 23d36'19.52"E,  0d54'51.94"N)
Lower Right (   18357.040,    -941.200) ( 23d36'19.52"E,  0d 1'51.74"S)
Center      (   15748.720,   13393.640) ( 23d31' 9.83"E,  0d26'30.10"N)
Band 1 Block=10032x1 Type=Float32, ColorInterp=Gray
  Min=-0.004 Max=0.048
  Minimum=-0.004, Maximum=0.048, Mean=0.004, StdDev=0.003
  NoData Value=-3.40282265508890445e+38
  Metadata:
    STATISTICS_MAXIMUM=0.047876562923193
    STATISTICS_MEAN=0.0043105538663359
    STATISTICS_MINIMUM=-0.0035738591104746
    STATISTICS_STDDEV=0.0025277971613299


texts: Overlaying layers with different projections on Android

texts:
Hello, I'm trying to show a map on Android. I'm using two map providers and am having trouble when converting geo point(lat/lon) to point in my viewport(x/y in pixels).

1) the first map provider used wgs84 -> EPSG:3395 projection
2) the second map provider used EPSG:4326(wgs84) -> EPSG:3395 projection

In the  first case I can place any lat/lon correctly.
In the second case I have a discrepancy. I'm not sure about 2 km or more.

Why do I get this discrepancy if the projection of the providers is the same?

Please help me understand this.

UPD: This is the java code

        double
        d = lon * Math.PI / 180,
        m = lat * Math.PI / 180,
        l = 6378137,
        k = 0.0818191908426,// somebody explain what this constant means?
        f = k * Math.sin(m);

        double
        h = Math.tan(Math.PI / 4 + m / 2),
        j = Math.pow(Math.tan(Math.PI / 4 + Math.asin(f) / 2), k),
        i = h / j;

        double[] res = {l * d, l * Math.log(i)};
        return res;

This code works for wgs84 -> EPSG:3395.  I know thsi is the wrong approach but I just changed constant k = 0.0818191908426 to my constant and it should work for some geographic region, but I think on the pole it will be wrong.


texts: Getting average elevation in floating point format around NumPy point coordinate?

texts:
Based on How do I get the pixel value of a GDAL raster under an OGR point without NumPy?, I wrote the following python function to find elevation height values of given coordinates in a DEM file:

from osgeo import gdal,ogr
import struct
import numpy

def single_point_elev(dem_file, points):

    src_filename = dem_file

    src_ds=gdal.Open(src_filename)
    gt=src_ds.GetGeoTransform()
    rb=src_ds.GetRasterBand(1)

    rows = points.shape[0]
    res = []
    i = 0
    while i < rows:
        mx,my=points[i][0], points[i][1]  #coord in map units

        #Convert from map to pixel coordinates.
        #Only works for geotransforms with no rotation.
        px = int((mx - gt[0]) / gt[1]) #x pixel
        py = int((my - gt[3]) / gt[5]) #y pixel

        intval=rb.ReadAsArray(px,py,1,1,buf_type=gdal.GDT_Float64)
        # print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
        res.append(intval[0])
        i = i + 1
    nres = numpy.array(res).reshape(len(res), 1)
    return nres

if __name__ == "__main__":

    nodes = numpy.genfromtxt('tc_outlet_cor.txt')
    print single_point_elev('hc.tif', nodes)

and the result was:

[[ 107.]
 [ 109.]
 [ 108.]
 [  26.]
 [  25.]
 [  27.]
 [ 247.]
 [ 239.]
 [ 249.]
 [ 281.]
 [ 283.]
 [ 281.]]

In the case above, input DEM file was a region containing streams, but as I checked the result, I found that some height values which belonged to downstream section were higher than that of upstream points, and besides, I got two points had the same elevation height value, as showed in the result above.

How could I get average elevation height value around given coordinate in floating point format?


texts: Attribute table is missing in raster after using extract by mask tool

texts: Every time I use the extract by mask tool, the rasters that are created don't have an attribute table (i.e it's greyed out). Does anyone know what is causing this and why these rasters essentially no longer have any values associated with them? Is there a way to get around this issue? I am using ArcGIS for Desktop.


texts: Specialized symbology packages for QGIS?

texts:
I would like to display a map of different biotopes in QGIS. In ArcMap, I can use predefined symbology for various topics by searching for it; e.g. for nature conservation issues.

Is there anything like that in QGIS, or are there predefined symbology packages that I could download?


texts: QGIS - error connection of database in RiverGIS

texts:
Currently, I am facing a problem to connect Postgres database to RiverGIS. It shows this error.
Couldn't connect to river database: rivergis@localhost
Please, check you database connection settings!

But, I checked the connection in the PostgresSQL, and it says the connection is successful.


texts: EarthEngine: average/reduce a raster at resolution of coarser raster?

texts:
I want to average a 30m raster at the resolution of a 1000m one. Specifically, I want to computer, for each Daymet raster cell, the percentage of cultivated land from the 30m CDL.

How can I do that in EarthEngine? My question is very similar to one asked on the EE list. The author was exploring 3 solutions:

A bilinear resample
Reproject
ReduceRegions, using a coarser raster that has been converted to features.

Someone gave an answer using reduceNeighborhood() to use a mean reducer, then reproject(), but it does not give quite the right result (see below, not right scale). How can I adjust that solution (I'm giving wrong scale parameters presumably?), or is there another solution that would be better?

Code

// DATA Import
var CDL = ee.ImageCollection("USDA/NASS/CDL")
var daymet = ee.ImageCollection("NASA/ORNL/DAYMET_V3").
  filter(ee.Filter.calendarRange(2015, 2015, "year"))

var aoi = ee.Geometry.Polygon([[[-100.501, 42.9819],[-100.501, 42.4045], [-98.39, 42.40457], [-98.39, 42.98192]]])

var CDL_2015 = ee.Image("USDA/NASS/CDL/2015").select("cultivated").clip(aoi)
var DYM_2015 = ee.Image(daymet.first()).select("tmin").clip(aoi)

// Operation
var image_frac=CDL_2015.eq(2).reduceNeighborhood({
  reducer: ee.Reducer.mean(),
  kernel: ee.Kernel.square(250,"meters"),
}).reproject(DYM_2015.projection().atScale(1000)).rename("cultivated")

//Visualization
Map.centerObject(aoi, 12)

Map.addLayer(DYM_2015.randomVisualizer(), {},  'DAYMET tmin')
Map.addLayer(CDL_2015.select("cultivated").eq(2), {min:0, max:1, opacity: 0.8, palette: ["beaed4","7fc97f"]}, "CDL coverage original")
Map.addLayer(image_frac.select("cultivated"), {min:0, max:1, opacity: 0.4}, "CDL coverage")


texts:
I have been trying to find an effective method (be it an equation or table filled with values) that would be useful for converting differences in longitude values into a distance in kilometers.

It was easy to find for Latitude since 1 degree is more/less equivalent to 111.19km, so any difference between 2 Latitude values would be as simple as:

distance = 111.19*(latitude1 - latitude2)

However, Longitude varies when going to a different Latitude line. At the highest and lowest latitude values (90 and -90 respectfully), each degree of Longitude is 0km in distance (obviously since all Longitude lines meet up and are right next to each other).

At the Equator, 1 degree of Longitude equates to 111.32 km (obtained from a random online source).

Is there an equation that could be used to calculate the distance between two longitude points? I have been searching everywhere but have yet to find a source. And if someone believes a table would be better than an equation, would you happen to have a link to such a table?


texts: Using WMS which is private but publicly available. Is it legal?

texts: I recently came across a useful WMS with ranges of various cellular technologies. It is used on a webmap of a private company but a simple look into network activity allows you to see the service address. There is no mention whether the underlying service is for public use or not but it is unprotected. Is it legal to use such a WMS for own purposes?


texts: Accessing @atlas_geometry in QgsExpression in QgsFeatureRequest?

texts:
It seems evaluation of @atlas_geometry is always 'null'

test code:

exp = QgsExpression('if(@atlas_geometry is null,1,0)') // always returns 1
rq = QgsFeatureRequest(exp)

@atlas_geometry within layer's expressions seems work fine though.

I have to do double spatial filtering across two layers so I come up with writing a predicate in PyQGIS for using in layer expression. Unfortunately it seems @atlas_geometry always returns null in my python code, despite it works fine in layer expression.

Can I access @atlas_geometry in my python code as if it's in layer's expression?


texts: GeoServer property interpolation datastore

texts:
I have a problem concerning our GeoServer and implementing a datastore.
Regarding the docs https://docs.geoserver.org/latest/en/user/data/app-schema/property-interpolation.html
it should be possible the use variables in the definition of a datastore.
I am trying to get a Oracle NG connection to a database using the web-interface.
When inserting normal values for the variables it is working.
When I shutdown the GeoServer and change the parameter database to
${example.database} in the datastore and starting it back up it is not working, I just get a null value when trying to add a layer from that datastore in the web-interface. When I look inside the file via web-interface there is still the value ${example.database}
My App-Schema Properties resides in
/data/tomcat/webapps/geodata/WEB-INF/classes/app-schema.properties

In the Tomcat startup I defined:
-Dapp-schema.properties=/data/tomcat/webapps/geodata/WEB-INF/classes/app-schema.properties

inside app-schema.properties: example.database = databasename
I checked permission for folders and the file itself.
Do I need to create the datastore manually and not the way I am doing it currently via web-interface and then inserting the parameter with an editor?
Is there a problem because my webapp is not named GeoServer?
Or is there something I'm missing like renaming the app-schema.properties or placing it directly in the datastore folder?
This is the datastore.xml:
<dataStore>
  <id>DataStoreInfoImpl--78867253:1790e655e57:-7fff</id>
  <name>test_db</name>
  <description>test_db</description>
  <type>Oracle NG</type>
  <enabled>true</enabled>
  <workspace>
    <id>WorkspaceInfoImpl--4db88a07:1773e4d4c5e:-7fb8</id>
  </workspace>
  <connectionParameters>
    <entry key="Evictor run periodicity">300</entry>
    <entry key="Estimated extends">true</entry>
    <entry key="Max open prepared statements">50</entry>
    <entry key="fetch size">1000</entry>
    <entry key="Expose primary keys">false</entry>
    <entry key="validate connections">false</entry>
    <entry key="Connection timeout">20</entry>
    <entry key="Batch insert size">1</entry>
    <entry key="Metadata bbox">false</entry>
    <entry key="database">${example.database}</entry>
    <entry key="port">1521</entry>
    <entry key="passwd">the_passwort</entry>
    <entry key="min connections">1</entry>
    <entry key="host">variable from etc/hosts which is working</entry>
    <entry key="Loose bbox">true</entry>
    <entry key="dbtype">oracle</entry>
    <entry key="namespace">http://xxx/geoserver/allgemein</entry>
    <entry key="max connections">10</entry>
    <entry key="Evictor tests per run">3</entry>
    <entry key="Test while idle">true</entry>
    <entry key="user">the_user</entry>
    <entry key="Max connection idle time">300</entry>
  </connectionParameters>
  <__default>false</__default>
  <dateCreated>2021-05-11 07:36:25.152 UTC</dateCreated>
  <dateModified>2021-05-11 07:39:24.759 UTC</dateModified>
</dataStore>

The error when trying to create a feature out of this datastore:
Caused by: java.lang.RuntimeException: Unable to obtain connection: Cannot create PoolableConnectionFactory (Listener refused the connection with the following error:
ORA-12505, TNS:listener does not currently know of SID given in connect descriptor


texts: Missing 'values' picklist in QGIS expression builder

texts:
Recently converted to QGIS 2.18.2 from Mapinfo under Windows10.
The Expression Builder in all the menus I've tried fails to display a popup window of values under the 'fields and values' selection when I click or right click on a selected field within the Expression builder window. Doesn't work for either 'top 10 unique values' or 'All unique values'.

However, if I right click a layer in the layers panel and choose 'filter' I can get pick lists for each field.
I've uninstalled QGIS, removed all the plugins, re-installed, still no joy.

I figure I'm missing something simple, but have worked through lots of tutorials. Am I missing something obvious?


texts: Combine two collections (MOD13Q1 and MYD13Q1) from Google Earth Engine

texts:
Need to combine two image NDVI collections (MOD13Q1 and MYD13Q1) in Google Earth Engine but don't know how to proceed.

My proposed workflow is to combine the two NDVI products from two collections to create a multi-temporal dataset. Then, from there will create a stacked image (image layer stacking) composed of NDVI images from the combined dataset, which will be later used for classification.

var region = ee.Geometry.Rectangle(xxlon, xxlat, xxlon, xxlat);
var roi = region;

// Load MODIS MOD13Q1 dataset.
var MOD13Q1_fullCollection = ee.ImageCollection('MODIS/MOD13Q1').select('NDVI');
var MOD13Q1_dateFilter = MOD13Q1_fullCollection.filterDate('2015-01-01','2017-12-31');
var MOD13Q1 = MOD13Q1_date.filterBounds(roi);
print (MOD13Q1);

// Load MODIS MOD13Q1 dataset.
var MYD13Q1_fullCollection = ee.ImageCollection('MODIS/MYD13Q1').select('NDVI');
var MYD13Q1_date = MYD13Q1_fullCollection.filterDate('2015-01-01','2017-12-31');
var MYD13Q1 = MYD13Q1_date.filterBounds(roi);
print (MYD13Q1);

var stackCollection = function(collection){
                      var first = ee.Image(collection.first()).select([]);
                      var appendbands = function(image, previous){
                      return ee.Image(previous).addBands(image.reproject('EPSG:4326',null,232));
                          };
                    return ee.Image(collection.iterate(appendbands,first));
    };

var stacked = stackCollection(supposed_to_be_the_combined_NDVI_from_two_collections);


texts:
I have a PostGIS database that gets updated about once every two hours. Since Geoserver datastore do not do any caching, I was looking for service does that and i found GeoWebCache REST API http://geowebcache.org/docs/current/rest/seed.html

from what i understand, its the only option for caching updated (new) data. am i right? and is there any detailed example on how to apply it?

I am using openlayers, with wms


texts: script for automated line drawing procedure for set of start and end points

texts:
A colleague conducted a field survey and recorded the start and end points of the transects with a GPS.  These points have been imported to ArcMap 9.3 as 2 point layers.  Both layers have a "transect number" field where each pair shares the same value, in order to match up start and end points.

I can manually "connect the dots," but there are about 1000 pairs of points, and I am hoping there's a quicker way.

Does anyone know of a simple script to run this procedure?


texts: Sharing QGIS Project with PostGIS connection

texts:
I am trying to share a QGIS project with multiple people which accesses a PostGIS database. Is there a relatively simple way to do this? How can I make sure that they are able to access the database?
Edit: If I share a QGIS project with a few users outside my network, will they be able to access the project if the project connects to a local PostGIS database on my computer, or will I somehow have to provide them access to the database? I am new to QGIS and could use some help. I have read some posts about editing the pg_hba.conf file, using GeoGig, or using AWS to host the database, but I am still confused as to which approach to take and why.


texts: Importing Excel point and line data into ArcGIS Desktop as shapefile?

texts:
I have two separate excel documents that I'd like to import into ArcMap, but I don't exactly know how to create a shapefile. The two files are as follows:

Point data. This spreadsheet has a long list of points (road intersections) with unique IDs and XY coordinates. This one shouldn't be too hard to import.
Line data. This spreadsheet has "From" and "To" columns with the unique IDs from the point data spreadsheet. It also has the XY coordinates of the points.

Is there any way to import the line data spreadsheet so that ArcMap draws a line between each of the points?

Another thing to consider is that many of these roads are not simply straight lines, so I don't know how to account for that. I'm using ArcGIS Desktop 10.1.


texts: Diminishing performance of gdal_polygonize in a python subprocess for GRASS to PostgreSQL

texts:
I have written a python script to polygonize GRASS rasters into a PostgreSQL database using the gdal utility gdal_polygonize:

sub.call('gdal_polygonize.py -8 gdb/quercus/{0}/cellhd/{1} -f "PostgreSQL\" PG:\"dbname=quercus\" {2}'.format(r,rcRast,pgTableName), shell=True)

The command works, however performance seems to drop off as the process proceeds. I am woefully ill-prepared to debug performance issues that may be related to the use of subprocess.

Are there any suggestions for a more robust use of subprocess and/or ways of benchmarking the performance?


texts: Creating a menu interface to update an sde feature class in ArcMap using arcpy?

texts: I work in electrical utilities and we have stacks of paper records that need to be updated into an existing SDE feature class attribute table. My problem is that someone in the city's GIS office would like me to create a simple interface where anyone with access can search a for a pole number and edit records in that row of the attribute table with a menu-like interface. He/she's reasoning is that employees with little GIS experience won't make mistakes with the data. Is it possible to construct something like this that would provide any benefits over starting an editing session and just keying values into the attribute table?


texts: How to configure base layers in Geonetwork

texts:
I am trying to get an instance of Geonetwork configured but I can't get base layers to display in the Map Viewer as shown below:

I think problem has to do with the fact that I installed Geoserver before Geonetwork so the default WMS base layers that Geonetwork uses where not set up.  Is there a script somewhere that will add these or do I need to do it manually?


texts:
I am beginner and I'm having troubles to do what I want in Earth Engine.
What I want: put the coordinates under the map (under the layer). Like this:
I already made the layer, but I don't know how to put the coordinates there. What I already did:
var Afa = ee.Geometry.Polygon(
[[[-47.405, -21.911],
  [-47.271, -21.912],
  [-47.250, -22.017],
  [-47.421, -22.018]]]);
Map.centerObject(Afa,12);
var proj = ee.Projection("EPSG:22523").atScale(4000)
var coords = ee.Image.pixelCoordinates(proj)
              .clip(Afa)
              .convolve(ee.Kernel.gaussian(5, 3));
Map.addLayer(coords.subtract(coords.round()).zeroCrossing().reduce('sum').selfMask(), {},'Coordenadas', false)

Any idea?


texts: I am using FME workbench 2018.1. I am working on a model where i have to apply lots of string operations on different feature classes and then merge them into one. I want to know if there is any console area where i can quickly test the function and see the output rather then running the whole FME again and again. I have been using inspector and logger but i am wondering if there is a console in FME.


texts: Is it possible to change distance in links table in AequilibraE plugin?

texts:
I'm in a project, which want to compute fastest path for container truck with some path are blocked. I am using QGIS with plugin AequilibraE (get data from OSM) to solve it, but I can't describe blocked path in data. After that, I think if I change the distance (in database) of blocked path to very large, the algorithms compute the fastest path never go through it. But, I found this in the docs of AequilibraE plugin :

3.4.2.3.1.3.1. Link distance   Link distance cannot be changed by the user, as it is automatically recalculated using the SpatiaLite
function GeodesicLength, which always returns distances in meters.
http://aequilibrae.com/python/latest/project_docs/network.html

So is there a way to change it at any cost, or is there a better way to describe blocked path?
P/s: AequilibraE using SQLite database.


texts: QGIS raster image leaving original colors intact

texts:
I'm using a program that generates a map using a set of colors (radio signal strength ranges, to be specific). It outputs the file as a jpg image with KML, so it displays perfectly well in Google Earth and in any photo viewer. But when I try to read it into QGIS (I've got 3.2 installed now), the raster processing functions seem to demand that they take over and replace the colors. Using palleted unique values, I was able to painstakingly regenerate something close to the original, but none of the options seem to just want to display the colors as the jpg file shows them.

Is there a trick I'm missing?


texts: Which ArcGIS Online basemap has the best performance in ArcMap?

texts:
My company plans to do demos of ArcMap Add-Ins for some of our clients and we want to make sure ArcMap is performing the best it can.

One of the major performance concerns is the speed of the basemaps found under
File-> Add Data-> Add Basemap.

Right now, we usually use the "Light Gray Canvas" basemap because of its simple appearance but recently I tried the more complex looking "National Geographic" basemap and it seemed to run better. I believe the basemaps are provided by different institutions so their IT infrastructure might have an impact on load times.

So I am wondering which basemap would have the best performance in ArcMap. I have already looked into things such as hardware acceleration, SSDs, and the cache settings.


texts: Removing whitespaces with layerlist and arcpy.Describe

texts:
This is an old problem with a new twist. I need to take map packages and extract all of the temp data into a geodatabase using feature class to feature class. I need a good work around for getting rid of spaces in baseNames before sending them through arcpy.Desribe(x).baseName there are simply too many files in too many temp locations to do this the hard way. Here is my naming script, before I throw my feature class to feature class function script:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
lyrs = arcpy.mapping.ListLayers(mxd)
for layer in lyrs:
    name = arcpy.Describe(layer).baseName
    name = name.replace("-", "_")
    #name = name.replace(" ", "_")
    layer.name = name
    print name
arcpy.RefreshTOC()

Obviously the name.replace(" ", "_") isn't going to work because it is after the arcpy.Describe function. Is there a way I can use arcpy to make a list of the layers and use that list with os to access each layer in a loop to change the filename? One problem with this solution is that half of the mpk files are in a geodatabase and half are lone shape files. Obviously those in a gdb are in correct formatting though.


texts:
I know there are a number of similar questions, but none of them have answers that allow me to use shapely in QGIS. In particular I am trying to install the Contour plugin, but this fails with the message "The plugin depends on some components missing on your system. You need to install the following Python module in order to enable it: shapely.geometry." The reason for this error is explained well here: How to install the Python module shapely for QGIS use?

These are the steps I have taken to solve the issue:

I first installed shapely via pip. This allowed me to import shapely in python 2.7, which is the default on my system, but installing the plug in still failed with the same error.

This question How to install shapely for python 2.6 (Mac) NOT 2.7? explained that qgis uses python 2.6, not 2.7, so I downloaded the tar file and did a  sudo python2.6 setup.py install. Now I can import shapely in python 2.6 too.

But the installation of the plugin still fails with the same warning! I've also tried import shapely from the python console within QGIS, and there is "no module named shapely." Is there something I've overlooked?

EDIT. Thanks to the comments from gene below, it turns out that QGIS on Lion does use python 2.7, not 2.6. But still it doesn't work. It turns out though that I have a few different versions of 2.7:

From the terminal:

$ python
Python 2.7.2 (default, Mar 8 2012, 15:29:59)
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.1.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>

From the QGIS python console:

>>> import sys
>>> print sys.version
2.7.1 (r271:86832, Jun 16 2011, 16:59:05)
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)]
>>>

From the terminal again:

$ python2.7
Python 2.7.3 (default, Oct 22 2012, 06:12:32)
[GCC 4.2.1 Compatible Apple Clang 3.1 (tags/Apple/clang-318.0.58)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>

So it seems that QGIS uses 2.7.1, Apple defaults to 2.7.2, and I also have 2.7.3. So I need to install shapely for 2.7.1?


texts: How to avoid the faint display of tracks when using the time slider of Google Earth?

texts:
After playing a track with the Time Slider, I often get a faint but clearly visible track remaining. Anyone know why and how to fix it please?

Here's an example. The GPX track has been played once and now I'm stuck with this faint remainder.

https://dl.dropboxusercontent.com/u/4019461/GE-TimeSliderFaintPuzzle.jpg


texts:
I´ve just tried to install "Profile tool" plugin on my QGIS version (2.14) but the following error message appeared:

"The plugin is broken. Python says: No module named PyQt.QtCore".

Does anybody any idea about how to solve this problem?


texts:
I'm currently using CDO to extract data from netCDF files containing satellite meteorological data. In the netCDF file points have a resolution of 0.5 degrees, so when data for a specific coordinate x,y is requested, I use cdo -remapbil,lon=x/lat=y to generate a temporary netCDF file from which to extract the data relative to x,y. Now as new data is being added the process is getting slower and slower.

In search of new ways to accomplish the task, I see that PostGIS can save raster data (in fact I've been able to import some data, converting netCDF to GeoTIFF then using raster2pgsql to import the data).

Now I'm stuck because I see that the raster is imported as pixels and I don't know how to formulate a query for a generic x,y coordinate that would return the same data as with CDO. How shall I tackle it?


texts: Seeking software for making digital elevation models from UAV aerial imagery

texts:
I have an unmanned aerial vehicle (UAV, also called a drone) and would like to use it to create a digital elevation model and orthomosaic over a smaller area (10 km2). My platform will be able to produce overlapping imagery with an altitude of around 120 meters (400 feet).
Several commercial products seem to provide complete solutions; eg. Pix4D and Dronemapper.
I am however primarily looking for solutions where I can have more control of the process and possibly also set up my own batching processing. Python and other open source tools are interesting.
What are my alternatives?


texts:
I processed one country with TPI, and I got a certain Landforms classification of terrain.
(Topographic Position Index (TPI) calculation as proposed by Guisan et al. (1999)).
When I processed only one municipality around that state with TPI, I got a different result, in fact more detailed. The calculation parameters were the same. This worries me because I thought to analyze the whole country, and then I would use the zonal statistics tools for each municipality individually and make conclusions based on the obtained forms of the field.
I am writing a scientific paper in this field, so I really care about real data.


texts: How to choose with a different way the OSM Boundaries

texts:
I need to get the borders of a country via Overpass API, but with the following query

relation
  ["boundary"="administrative"]
  ["admin_level"="2"]
  ["name:en"="some country for example Papua New Guinea"];
(._;>;);
out body;

I get the marine borders also, and not only the borders of the terrain I want, any solution for that?


texts: How to get 2010 city demographic data by census block group

texts:
I am doing machine learning research on US cities (New York City, Los Angeles, San Francisco, Chicago, Seattle). I try to predict some demographic data based on features in each census block group.

I am looking for 2010 household median income data by census block group in five cities above. I hope the data can contain the census block group location (lat/long or coordinates) and its median household income. I have looked through some database websites such as FactFinder, NYCOpenData, US Census Bureau, and I also looked through similar questions....but I haven't found what I want.

Any idea where I can get the desired data? Or is it possible to get such desired data?


texts: How to Set up Data Driven Pages using arcpy?

texts: I've read a lot about the topic data driven pages. I'm writing a python script for ArcGIS and  I've defined my stripmap using the arcpy.StripMapIndexFeatures_cartography command and now I have to set up Data Driven Pages with all the parameters it needs. My only problem is that I cannot understand how to do this. Is there a specific command in the arcpy library or a specific workaround to do this?


texts: python script in modelbuilder: get field value and return boolean expressions

texts:
I have problems integrating a python script in a modelbuilder tool. I'm a absolute beginner with python. Maybe there's also an easier way to solve this problem in modelbuilder without using a python code, but I haven't found it yet

so here's my problem:

my modelbuilder tool's flow of control needs to evaluate a situation and decide on one of two courses of action: If the value in the field "Building" in my feature class is 1, it should follow one path, if the value is greater than one it should follow another path.

Here's my python code that I tried so far, but didn't work.

# Load the arcpy module
import sys, os, arcpy

# Get the input feature class from the model
InputFC = arcpy.GetParameterAsText(0)
field = "Building"
cursor = arcpy.SearchCursor(InputFC)

# check if one or more than one building footprint exist
for row in cursor:
    Building =(row.getValue(field))
    if Building == 1:
        arcpy.SetParameterAsText(1, "True")
        arcpy.SetParameterAsText(1, "False")
    else:
        arcpy.SetParameterAsText(1, "False")
        arcpy.SetParameterAsText(2, "True")

any help on this problem is much appreciated

regards nicolas


texts:
I'm having difficulty with a PostGIS query. I am trying to query my database (world_countries) using a bounding box drawn on a google web map. The objective of this query is to find any countries that lie within a box drawn by the user, and to return them in GeoJSON format. The query is as follows:

SELECT ST_AsGeoJSON(the_geom)
 FROM
 "world_countries"
WHERE
((the_geom && 'POLYGON((
-1572563.9181513 7593219.6338481,
-1709539.0728192 6507202.3361235,
-413167.07328311 6477850.5172661,
-271299.948805577549191.905562,
-1572563.9181513 7593219.6338481))'));

and the error message is:

ERROR: Operation on two geometries with different SRIDs
SQL state: XX000

The SRID of this data is 900913 (srid(the_geom) = 900913). When I look at this data in ArcGIS and pan around, the values from the box above match up with ArcGIS.

Is there a reason why a polygon geometry requires a specific SRID? I assumed that my values were within the range.

A quick note about my projection. I originally had a map of countries of the world in geometric form using WGS_84. I then used ArcGIS to re-project this into 'WGS 1984 Web Mercator.prj' based on this thread.

As I am using the following version of PostGIS "POSTGIS="1.3.6" GEOS="3.1.1-CAPI-1.6.0" PROJ="Rel. 4.7.1, 23 September 2009" USE_STATS" I ran the following SQL to include the google projection info:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) values (900913 ,'EPSG',900913,'GEOGCS["WGS 84", DATUM["World Geodetic System
1984", SPHEROID["WGS 84", 6378137.0, 298.257223563,AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], NIT["degree",0.017453292519943295], AXIS["Longitude", EAST], AXIS["Latitude", NORTH],AUTHORITY["EPSG","4326"]], PROJECTION["Mercator_1SP"],PARAMETER["semi_minor", 6378137.0],
PARAMETER["latitude_of_origin",0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor",1.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0],UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH],AUTHORITY["EPSG","900913"]] |','+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
+nadgrids=@null +no_defs');


texts: Connect to local file or PostGIS layer in field calculator QGIS

texts:
I'm trying to make query from field calculator so that it would be possible to connect to PostGIS layer (or local file on computer) and select data according to condition so that use this selection afterwards for layer in project.

Is it possible?

I've found great functions in Reffunctions plugin - e.g. dbvalue, that retrieve first targetField value from targetLayer when keyField is equal to conditionValue;

dbquery, that retrieve first targetField value from targetLayer when whereClause is true (for file in project) - this function seems to not work

But all these functions require a layer in a project, altough I wonder to connect to PostGIS layer or local file in a query. I suppose I will need to create custom function in function builder, but I don't know how turn to local or PostGIS layer from QGIS via python.

This post and cookbook gives information about adding PostGIS layer in QGIS project, but still unsuccessfully - What's wrong? It gives me an error

"NameError: name 'self' is not defined"

other variant

also don't add a PostGIS layer

I tried single quotes instead of double, but it doesn't help. Local file adds successfully after replacing single delimiter to double in path to file.

Any ideas?)


texts: Calculate distance of point along a line using NetTopologySuite

texts:
I am looking to calculate the distance along a line of a point using NetTopologySuite.

Indications are that the distanceOp may have something to do with it, given that I can a point on a line closest to another geometry.  However how to determine the distance along the line of that point?  I can certainly determine what that distance is myself through some code - but I would have through NTS would already have this covered off.

Anyone have any pointers?

Thanks,

Matt


texts:
I have a script that has several lines and include 3 arcpy functions such as copy features, make xy layer event, and apply symbology. Right below those tools I am trying to use the arcpy.mapping to create a mxd file with the new layer and export as JPEG. However, when I do that it seems the script doesn't follow its order and the output (the mxd file and the jpeg) don't have the new layer, which was created and analyzed during the script.
*obs: I am trying with the script tool and it doesn't work, when I try just putting the script in the ArcMap python window, it works perfectly fine.

I am having problem on this part:

import csv
import os
import string
import arcpy
from arcpy import env
env.workspace = "E:/Alex"
#variables
inp = arcpy.GetparameterAsText(0)
output = arcpy.GetparameterAsText(1)
#directories
outfile = "pathname"
jpeg = "pathname"
finalmap = "pathname"
    #csv reader
with open(inp, "r") as csvData:
    reader = csv.reader(csvData, delimiter=',', quotechar='|')
#csv writer
    with open(outfile, "wb") as csvWriter:
        writer = csv.writer(csvWriter, delimiter=',', quotechar='|')
try:
    arcpy.env.overwriteOutput = True
    outlayer = "Event_Layer"
    dataset = "examplelayer"
#Getting coordinate system from another layer
    spatialreference = "GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 11258999068426.2;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision"

    arcpy.MakeXYEventLayer_management(outfile,"lon","lat",outlayer,spatialreference,"")
 except Exception as err:
    print (err.args[0])

#Copying Features
arcpy.CopyFeatures_management(outlayer,output)
#Map Layout
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "")[0]
for lyr in arcpy.mapping.ListLayers(mxd):
    arcpy.mapping.AddLayer(df, lyr, "TOP")
    arcpy.ApplySymbologyFromLayer_management(lyr, dataset)
arcpy.mapping.ExportToJPEG(mxd, jpeg, "PAGE_LAYOUT")`
mxd.saveACopy(finalmap)
del mxd, outlayer


texts: Converting hours to acceptable FME format?

texts:
I have an Excel document with a 'time' column. I need to get this column into an appropriate format for FME to convert the values to minutes.
The problem is that the cells are in different formats (h:mm, [h]:mm:ss).

An example of the data is:

1:44
0:34
45:44
223:44

I.e 223 hours and 44 minutes needs to read 13424 minutes.

Any advice on transforming these values?


texts: Fill not filling all sinks using ArcGIS Spatial Analyst and ModelBuilder?

texts:
I am trying to make a model where I get a shapefile of conical hills in ArcMap 10.5.1. First I invert the DEM to get depressions. Then I used the fill tool and I set z-limits for each meter from 10 to 150. The problem is that the depressions on the higher altitude of my research area aren't recognised til z-limit = 150, not each depression seperatly but all of them together in one big depression with multiple sinks.

Is there a better way of detecting conical hills or is there a possible solution for my problem?

Here is a picture of my inverted DEM with some polygons. As you can see many are missing. I did not include the polygon produced at z-limit 150 so that the problem is visible.

And a picture of my model.


texts:
I am new to using OGC services and am attempting to extract a large tiff using the following URL request:

http://gis.dnr.alaska.gov/terrapixel/cubeserv/ortho?DATASTORE=SDMI_ORTHO_RGB&SERVICE=WCS&REQUEST=GetCoverage&
ACCEPTVERSIONS=2.0.1&coverageId=image_set&
SUBSET=E(101688,206499)&SUBSET=N(1005498,1100573)

At this extent, the resulting tiff would be ~13G, exceeding the server's max of 2G per request. However, when I try incrementing the extent to be withing the download bounds, I'm getting a 504 gateway timeout. I'd imagine that's probably still due to requesting a large amount of data, but I'd rather not go down the route of looping through 100 URL requests.

Obviously, I'd rather download the entire extent in one go, to avoid merging the rasters together, so ideally I can just add a new parameter that compresses the image. Like I said, I'm new to using WCS and not sure of the scope of capabilities/best practices so any advice would help a ton.


texts:
I have a database with geometries encoded in a field in the following format:

'1  L 410012.85,416177.69; +30.16,4; 12.46,0; 13.35,-0.67'

Where L specifies that it is a Polyline, then 410012.85,416177.69 is the starting point. The rest of vertices are calculated by adding or subtracting X and Y. So for example:

Vertex 1 is starting point (X=410012.85 Y=416177.69)
Vertex number 2 is X=410012.85+30.16 Y=416177.69+4
Vertex 3 is Vertex2(X)+12.46 Vertex2(Y)+0
Vertex 4 is Vertex3(X)+13.35 Vertex3(Y)-0.67

Has anyone seen something like this before?

Generating a standard format by splitting the text in multiple parts is something that probably can be done. However, I'd be more interested in knowing if it is actually a known format in GIS (does it have a name?) and if already exists a way of translating it into a standard format.


texts: How to record tracks finer than 10-30cm resolution (possibly using an iPhone)?

texts:
I am trying to record tracks at a much finer scale than is possibly with GPS but hopefully using something simple and cheap, ideally an iPhone. The resolution for track data I am after is around 10-20cm.

The location of the track, and even the scale is not as important as the fine details of the shape of the track.

What I have in mind as the ideal solution is an iPhone app that uses the camera in conjunction with the other sensors to record a track in a similar way to an optical mouse detecting movement over a desk's surface.

So you would hold the phone with the camera facing down and walk the track. The app would trace the moving surface and infer where it is going, using the gyro and accelerometer to help correct for tilting the camera. It would also use the gyro and to a lesser extent the compass to help correct and detect the turns in the track. The height of the phone off the ground I don't think could be sensed or inferred but you could just enter that and hold it carefully at a consistent height. The led could be on to aid detection through shadows.

If you happen to be walking in a track that finished where it started you could mark that in some way to close the loop and it would try and correct it if the ends didn't meet up, and if you walked a long way, perhaps > 50m it could also use gps data to calibrate the scale of the track, but probably not really contribute much to the details of the tracks geometry.

So does something like this even exist in any form iphone or otherwise? Is there anything which is simpler eg it just uses the gyro/compass and a pedometer style measurement method? Is there some completely different cheap simple method?


texts: Making Map Book with facing pages numbered near outside edge?

texts:
I am creating a map book but do not want the page numbers for each map page to be on the same side (since its printed double sided), I need the page numbers to be placed on opposite sides.

I am using data driven pages.

How do I do this?


texts:
I need to create image collection of NDWI for a polygon I have from 1983-2019.
For that, I want to use different Landsat image collection.
I saw that there is made a product in Google Earth Engine that gives the NDWI values, but I face few problems, the first one is that I want to represent and download the whole image collection, and the second is that whenever I try to represent even one image I don't see anything. I tried also to represent only one but I still see nothing.

This is the code I wrote.

// create image collections for the years 1983-2019

//Landsat 4
var dataset = ee.ImageCollection('LANDSAT/LT04/C01/T1_8DAY_NDWI')
                  .filterDate('1983-01-01', '1984-12-31')
                  .filterBounds(table);
print('dataset: ', dataset);

// Get the number of images Landsat4
var count = dataset.size();
print('Count: ',count);

// Sort by a cloud cover property, get the least cloudy image- why does it rank the images? do they still have clouds?
var image = ee.Image(dataset.sort('CLOUD_COVER').first());
print('Least cloudy image: ', image);

var colorized = dataset.select('NDWI');
var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: ['0000ff', '00ffff', 'ffff00', 'ff0000', 'ffffff'],
};
Map.centerObject(table);
Map.addLayer(colorized, colorizedVis, 'Colorized');
Map.addLayer(image,colorizedVis,'Hi');


texts:
In ArcGIS the function called "Auto Registration" and it can automatically georeference raster

Is there a way in QGIS to do the same?
I seek in the transformation settings and in the setting georeferencer but found nothing.


texts: Determine points where a LineString meets the boundaries of a Polygon using Turf.js

texts:
Using Turf, I'm trying to determine the points at which a line (LineString) meets the boundaries of a polygon.

Turf has the ability to do point-in-polygon/point-on-line and polygon intersections etc, but I can't figure out a way to get Turf to give the answer to the question "what are the points at which this LineString meets the boundaries of this Polygon".

The only options I've come up with so far is to enumerate the points of the polygon as line-segments and then use a line intersection algorithm.  I just hoped there was a more "Turf" way of doing it.

PS - I already asked this over on stackoverflow but someone suggested I try here on gis.stackexchange.com


texts: Roads overwrite railways in deep tunnels

texts:
I'm trying to set up a tileserver to deliver a custom map of London's railway system using CartoCSS and Mapnik. London has overground lines, underground lines, a light railway and a tram system.
I have split the railways into a separate layer from the roads and this railway layer is the last layer in project.mml. If I comment out the road layers, the railways are drawn correctly. If I comment out the railways, the roads are drawn correctly.
However, when I draw both the road and railway layers, the bits that run in tunnels are displayed correctly only where they run underneath parks and water; where the tunnels are under roads, the roads obscure the railways in the tunnels.
While there are aspect of attachments, the layernotnull field and z_order that I don't fully understand and might not have fully experimented with yet, my understanding is that each layer in the project file is drawn, top to bottom in the file (bottom to top on the map canvas) and so I don't understand how the roads always come out on top.
Am I missing something obvious?
Edit:

Here, the purple line is an underground railway, being drawn in the last layer in the project.mml file. It appear correctly on top of buildings and parks but roads are drawn on top of it - it is this bit I don't understand.


texts: I am looking for a way to add a picture to my map using pyqgis. Are there any examples of how to do this. I know it can be done using the classQgsComposerPicture. I am not much of a programmer so I cannot decode the api documentation. Any links where I can find such examples or documentation with options.


texts: Trouble updating attributes of a layer based on another layer's attributes if features intersect in python, QGIS

texts:
I am trying to assign attributes from features in one layer that intersect features in another layer. When I run my code, I do not get an errors but it does not update the attribute field as I am expecting. I wrote my code based on the example in this post https://stackoverflow.com/questions/23506358/pyqgis-for-spatial-join. It appears as though the program does not see that any features intersect, but when i compare them on the map canvas, they do. I'm not sure what's wrong. Here is my code:

GPSDegLyr.dataProvider().addAttributes([QgsField('UniqueID',QVariant.Int)])
GPSDegLyr.updateFields()

idxGPS = GPSDegLyr.fieldNameIndex('UniqueID')
idxBuff = dsslvExplodeLyrWGS84.fieldNameIndex('UniqueID')

GPSDegLyr.startEditing()

for b in dsslvExplodeLyrWGS84.getFeatures():
    for g in GPSDegLyr.getFeatures():
        if g.geometry().intersects(b.geometry()) == True:
            val = b.attributes()[idxBuff]
            print val
            g.setAttribute(idxGPS,val)
            GPSDegLyr.updateFeature(g)

    GPSDegLyr.commitChanges()

QgsMapLayerRegistry.instance().addMapLayer(GPSDegLyr)

UPDATE: I solved the intersecting problem - now the programs sees the intersections and prints the values. But still doesn't write the values out to the attrribute table.


texts: Merge Shapefiles and Add File Name to Attribute Table

texts: I am trying to merge lots of point shapefiles together to form 1 large shapefile. What i want to be able to do is get the name of each of the point files to appear in the attribute table either during the merge or carrying out a post processing task. How can i do this in ArcMap 10.3? I also have no experience with python code etc.


texts: Getting started with QGIS plugin development?

texts:
I wish to develop a simple plugin for QGIS, which will open some data, and perform some specific processing.

I tried searching on the QGIS wiki, but only found this link :

http://www.qgis.org/wiki/Writing_C++_Plugins

which seems to be obsolete. There were also a few blog posts from 2006.

I also saw something about Python plugins, but I am not sure whether I need to work with the C++ plugin or a Python plugin.

Where can I find documentation, or tutorials, that would help me get started?


texts: Querying polygon data imported from Google BigQuery GIS within R

texts:
I have a geography data type column that each row is a polygons (representing South Africa suburbs) which I have to extract from a BigQuery database into R and from there plot these in an R Shiny app (or use for other analysis).

Unfortunately, I keep running into this error:

  Downloading 77,850 rows in 8 pages. Error in
  bq_parse_files(schema_path, page_paths, n = page_info$n_rows,  :
  Unknown type GEOGRAPHY

I have tried various options such as:

query <- "SELECT *, ST_ASTEXT(ST_GEOFROMTEXT(wkt)) AS Geometry FROM od_matrix.cartographic_partitions"
tb <- bq_project_query("PROJECT_HARM", query)
data <- bq_table_download(tb) %>%
  sf::st_as_sf(wkt = "Geometry")

and

DBI::dbGetQuery(conn, "SELECT *, ST_ASTEXT(wkt) AS Geometry FROM `PROJECT_HARM-prod.geography_data` LIMIT 10")

I keep having the same errors.

Here are some of the sources, I have tried:

https://github.com/r-dbi/bigrquery/pull/376
https://github.com/r-dbi/bigrquery/issues/313
Manipulating polygon data imported from Google BigQuery GIS within R


texts:
I need to create a cloud mask using the pixel_qa band, this is the code that I have so far:

var city = ee.Geometry.Point(9.15, 45.18);
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

// Fmask classification values

var mosaic = landsat8
    .filterBounds(city)
    .filterDate('2017-08-01', '2017-09-01')
    .mosaic();

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

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

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

I tried to create a mask using the Map.addLayer but the result is only a white mask.
What I need is to create a mask in which clouds are displayed with a certain color (ex. white) and all the other pixel type with another color (ex. black).


texts: Using Select by Attribute with multiple layers in ArcGIS Pro?

texts:
I'm trying to find "Signal" in a field named "Maintained by" but there are quite a few layers.

I do not want to go into every single layer and change "Signal" to "Private."

Does anyone know a quick shortcut?


texts: Orthorectification with gdalwarp - projection issues

texts:
I'm working on gdal-based orthorectification workflow/procedure for WorldView2 images.  For my particular application, WorldView2 images will always cover somewhere in CONUS.  Here is a rough outline of this workflow:

Get WorldView2 imagery for AOI (area of interest)
Get 1/9 arc second (~10m) elevation data for AOI from USGS National Map
Convert elevation data to ellipsoidal heights (from NADV88 to WGS84).
Run gdalwarp on WorldView2 imagery for aoi with -rpc flag set, and -to (transformer option) flag set with RPC_DEM=dem_ellipsodial.tif (dem created in previous step)

Step 3 above is achieved using gdalwarp as shown below:
gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=g2012a_conus.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" dem.tif dem_ellipsoid.tif

I start running into issues at Step 4.  I initially tried calling gdalwarp like this:
gdalwarp -rpc -to RPC_DEM=dem_ellipsoid.tif wv2.TIF wv2_rectified.TIF

This usage of gdalwarp fails with the error:

Creating output file that is -2147483648P x -2147483648L. ERROR 1:
Attempt to create -2147483648x-2147483648 dataset is illegal,sizes
must be larger than zero.

Some further investigation of this error suggests to specify -s_srs and -t_srs flags.  In the particular area of I have been working in, the WorldView imagery and elevation data are in UTM Zone 17N.  So:

gdalwarp -s_srs EPSG:32617 -t_srs EPSG:32617 -rpc -to
RPC_DEM=dem_ellipsoid.tif wv2.TIF wv2_rectified.TIF

The above usage of gdalwarp completes without issue.  However, when I add the wv2_rectified.TIF into QGIS - it does not display where it should.  When I look at QGIS information of the wv2_rectified.TIF layer - I see:

Notice the extent coordinates.  It seems problematic to me that these coordinates are still in lat/lng even though the CRS shows as UTM Zone 17N.  What I am doing wrong here?
FWIW...  For the time being, I do have ESRI products available.  I have successfully used the "Create Ortho Corrected Raster Dataset" tool with the wv2.tif and dem_ellipsoidal.tif I am using in the gdalwarp command.  The output of this tool is good - so feel good about my conversion of dem data to ellipsoidal heights from NADV88.  I dont have access to ESRI for much longer, hence trying to figure out a gdal based alternative.


texts: Change variable for CartoDB torque heatmap

texts:
I've got a geocoded data set with campaign donations. I want to make a torque heatmap that reflects the size of the donations -- the point being that while there are more donations from within the state, most of the money is coming from outside the state.

I can create a torque heatmap that reflects the number of donations over time, but can't figure out how to have the data points reflect donation size rather than just the number of donations from that location. Any suggestions?


texts: gdal_translate with srcwin and center reference point

texts:
i want to use subset a large raster image using gdal_translate with srcwin parameter
but i want also to give some reference point to be at the center of the output file .

so if the output file size is 2000X2000 pixel   for exemple, i want that my given point(lat/lon) to be at the center

How can this be achieved?


texts: TypeError: can't multiply sequence by non-int of type 'numpy.float64'

texts:
I'm trying to do a regression, but it does not allow me to multiply my values within the exponential function.

import matplotlib.pyplot as plt
import scipy
import pandas as pd
import numpy as np
from scipy import optimize
import math

data = pd.read_csv("/Users/Bruno/PycharmProjects/Benzoic Acid/BA Basic
Data/Rot-Disc-016.xlsx-Table 1.csv")
len_d = len(data.iloc[2:,0])

t = []
Benz = []
DC = []
dC = []

for i in range(len_d):
    t.append(float(data.iloc[2+i,0]))
    Benz.append(float(data.iloc[2+i,4]))
    DC.append( float(data.iloc[1+len_d,4])-float(data.iloc[2+i,4]) )
    dC.append(float(data.iloc[2+i,5]))

def obj(beta, x, y):
    beta_0 = beta[0] #This will allow the elements to be converted into beta vectors
    beta_1 = beta[1]
    yhat = beta_0 * np.exp( x * beta_1 ) #This is the choosen model equation
    resids = yhat - y #Will calculate residues
    SSR = np.sum(resids**2) #Will produce the sum of the square residues(SSR)
    return SSR

result = scipy.optimize.minimize(obj, x0=[1,1], args=(t, DC )) #This
optimizes the intended function
beta_opt = result.x


texts: How to query all history of specific user's activity?

texts:
I would like to query all of the existing meta-data on a specific user using Overpass QL.  I am new to the format, but I know that this will return all the nodes in a box by a specific user:

node
    (uid:123456)
    (lat,long,lat,long);
out meta;

And I could run the same command for rel's and way's in order to get a complete data-set.  But this will not, as far as I can tell, show me the history of the user's activity, in that I won't be able to see the history of edits and changes, deleted items, etc.: only the stuff in the most recent changesets/versions.

I feel like there must be some command along these lines to accomplish what I want:

[diff:"YYYY-MM-DDTHH:MM:SSZ","YYYY-MM-DDTHH:MM:SSZ"];
(
  (uid:123456);
);
out body;
>;
out meta;

But as you can tell by running that here, it doesn't work.  Also, if it did work, it wouldn't really do what I want, since it would only tell me the difference between one date and the other, as opposed to a list of all the edits of any sort a person has done/made.  So the question is: how can I query all of the meta-data on the activities of a specific user in Open Street Maps?  (Doing so within a specific date-range and/or bounding box is a plus but not necessary, if your answer takes neither into account I can implement both features on my own after the fact).

EDIT -----
I was able to get results similar to what I'd like with this:

[out:csv(::type,::id,"name",::lat,::lon,::version,::timestamp,::changeset)]
[diff:"2012-08-14T15:00:00Z","2014-08-21T15:00:00Z"];
node
    (uid:1378531)
    (46.32701220,11.63249080,46.58063280,12.44677780)
out meta qt;
rel
    (uid:1378531)
    (46.32701220,11.63249080,46.58063280,12.44677780)
out meta qt;
way
    (uid:1378531)
    (46.32701220,11.63249080,46.58063280,12.44677780)
out meta qt;
 >;
out meta qt;

But as of yet I am not totally sure if these results have all of the changes, or just the initial version and the latter version of any edited node.  Moreover, when I try running it with larger date ranges, the server gives me an empty package.  I am not sure how I'm supposed to query when I want something that necessitates more memory. Should I break it into multiple date ranges if the range is more than, say, a month, and then combine the results on my local machine after the fact? (It's not hard to combine CSV files with bash so this is definitely an option...)

Per a comment, I should mention that I do not have the flexibility to download gigabytes of data and parse through them.  I need to be able to download only the relevant data.


texts: Optimized smoothing of stream lines generated from Arc Hydro?

texts:
Has anyone found an effective ArcGIS 10.3.1 workflow for smoothing streamlines generated from ArcHydro...or hydrology in general?

(Fill -> Flow Dir. -> Flow Acc. -> Stream Def. -> Stream Seg.-> Drainage Line Processing) outputs a polyline that requires manual QA/QC to determine if lines are real and flowing appropriately. Lines look pretty good at 1:1k.

The influence of the source 1 meter DEM becomes more apparent when zooming in.

ArcGIS smooth tool (Advanced Editing toolbar) allows the user to specify a maximum allowable offset. Although you can enter a value between 0 and 1, there does not appear to be any difference between 0.1, 0.2, 0.3, etc. and 1.

ArcToolbox -> Cartography Tools -> Generalization -> Smooth Line tool grants the user the ability to leverage the PAEK or BEZIER INTERPOATION algorithm. PAEK is the default within the editing environment. BEZIER appears better at capturing the far bends or oxbows in a meandering watercourse.

Additionally, Esri help documentation highlights the difference between shapefiles and geodatabase feature classes when smoothing. Only feature classes are capable of true Bezier curves in between adjacent vertices. Shapefiles approximate a curve with multiple straight line segments.

At the moment, I am deriving preliminary stream polylines from a 1 meter bare earth DEM (partially hydro-enforced) with the ArcHydro workflow above. I'm then comparing the output shapes to ancillary surface water data and editing the raster-derived lines manually to add or subtract vertices, so they line up with 6" resolution aerial imagery.

To avoid the jagged line appearance at large scales, my best guess is to:

Make sure all line editing is taking place on a GDB feature class.
Utilize Cartography Tools -> Smooth Line tool -> BEZIER_INTERPOLATION
Ensure close match to aerial imagery by editing the pre-smoothed feature class.

Any other tips for increasing accuracy/efficiency?

Is it possible to make the Smooth tool in the Advanced Editing Toolbar utilize BEZIER INTERPOLATION?


texts:
I am seeking to calculate the variable MASSSTAB through rounding the variable LAENGE and doing some basic math.

My codeblock currently looks like this:

x = long(round(!LAENGE!, -1))

def ifBlock(m, x, a= 50, b= 30):
  if x > 10 :
    m= ( x * b )
  else:
    m= ( x * a )
  return m

My Codeblock is:

enter code here

There is already a syntax error (parsing error) in the first line. Do I have to define "x" in another way?


texts: worldclim: "Current conditions (interpolations of observed data, representative of 1960-1990)"

texts:
I am interested in work with climate data since 2012, I found it in  http://www.worldclim.org/version1, but there is something draws my attention:

There it says: "Current conditions (interpolations of observed data, representative of 1960-1990)"

Does it mean the current data is something like a "temporal projection" based in data of 1960-1990?

or

Should I understand it like an interpolation of weather stations existing since 1960 with not including stations existing since 1990?


texts: Is it ok to group by geometry in PostGIS?

texts:
Is it OK to group by geometry?

We do a lot of counting points by polygon geometries, which involves intersecting the data first in order to count instances of students. by school boundary, block group, etc:

SELECT
  n.nbhd_id
  , count(*) AS count_burglaries
FROM denver.crime AS c
  JOIN denver.neighborhoods AS n
    ON ST_Intersects(c.geom, n.geom)

WHERE c.offense_ty LIKE 'BURG%'
GROUP BY n.nbhd_id

Second step is to wrap in a sub-query in order to join the geometry from the polygon table back to the count query:

SELECT

count.*
, n.nbhd_name
, n.geom
FROM (
SELECT
n.nbhd_id
, count(*) AS count_burglaries
FROM denver.crime AS c
JOIN denver.neighborhoods AS n
ON ST_Intersects(c.geom, n.geom)

WHERE c.offense_ty LIKE 'BURG%'
GROUP BY n.nbhd_id

) AS count

JOIN denver.neighborhoods AS n
ON count.nbhd_id = n.nbhd_id

But it seems you can also use geometry in GROUP BY:

SELECT
n.nbhd_id
, n.nbhd_name
, n.geom

, count(*) AS count_burglaries
FROM denver.crime AS c
JOIN denver.neighborhoods AS n
ON ST_Intersects(c.geom, n.geom)

WHERE c.offense_ty LIKE 'BURG%'
GROUP BY n.nbhd_id, nbhd_name, n.geom

Is including geometry in the GROUP BY a correct approach?


texts:
I have some GIS project to do, in this project I need to place something like 20,000 trees on a map.  So I did a web app with .net that includes map that based on OpenLayers with Bing base map.

My problem is that I need to get more closer to the base map, because the maximum scale that I can get is 50 ft, but I need to be much closer to place the tree in the right place.

Is there any way to make option to do more zoom? More then 50 ft?

I tried to change the zoom levels or the resolution, but nothing worked for me until now.  Maybe I need to use other base map then Bing?  Or maybe I cannot get closer on raster map?


texts: Spatial join in ArcMap vs. ArcPy (how to make them the same)

texts:
I am trying to replicate the join tool that is available in ArcMap by right-clicking on a layer, in an ArcPy script. Specifically, all I want is an attribute field of how many points are inside each polygon. Works perfectly inside ArcMap, but I must have it inside an ArcPy script.

My polygon file has 90 polygons, and my points file has 1,070 points. When I do the join inside ArcMap (screenshot below), my output is a polygon SHP with 90 polygons, and only a "Count_" field is added. (This is perfect.)

When I do it inside ArcPy, my output is a polygon file with 1,070 polygons, a bunch of extra fields, and no sum-count field. I've tried a few different things in the ArcPy script (i.e., "JOIN_ONE_TO_MANY", "CONTAINS", etc.), but no luck. What am I missing?

NOTE: I don't have the Arc license that allows for "summarize_within"

arcpy.SpatialJoin_analysis(Polygons.shp, Points.shp, outFeatureClass, join_operation = "JOIN_ONE_TO_ONE", join_type = "KEEP_COMMON", field_mapping = "", match_option = "INTERSECT", search_radius = "", distance_field_name = "")