texts: Truncating field names to 10 characters in ModelBuilder

texts:
I have a number of different animal GPS location points, and a bunch of different habitat characteristics. I have built a simple model that uses the Near tool to calculate the distance from each point, create a new field using Add Field, naming it after the habitat characteristic shapefile name (using Parse Path and inline variable substitution), and then populates that added field with the NEAR_DIST values. It then continues on for all of my habitat characteristics.
I built the model a while ago and have since changed the names of the habitat characteristics (and there are a bunch of them) and unfortunately in the renaming them they are all now above the 10 character limit for a field name. What I was wondering is if there was any way to truncate the habitat characteristic shapefile title to 10 characters so that the add field tool doesn't error out.
Examples of habitat names are - 2001 Upland, 2001 Harvest, 2001 Deciduous
Essentially I need to remove the year and the space and they will all be under 10 characters.
The tool I created runs through but doesn't successfully "Add Field" because it is trying to name it after the "Near Feature" (the habitats whose names are too long) and so it just skips it.


texts: QGIS 2.16 - trying to export PDF but minidump!

texts:
It just keeps minidumping.

Working on a big project and it just keeps mini-dumping, have even added more RAM to the computer to try make it easier.
So far today I've reinstalled 2.16 (64bit), deleted all the errors in the users temp file and tried removing 'QGIS2' from the file path. I recreated the project (using the layers I've already made) and have got so far as trying to export the PDF from print composer before POW 'minidump written to C: \Users\Belinda\AppData\Temp\qgis-20160824-121716-5792-6496-8545b3b.dmp'

How do I get around this? What else can I try?


texts: Stop QGIS on Windows from recentering on single-click

texts: I'm using QGIS on Windows 10. In Windows, the normal UI is that when you want to bring an application forward, you can click on it. And that tends to select what you clicked on, but single-click select usually doesn't do much; double-click executes. In QGIS, most of the screen area is taken up by the map. My problem is that single-clicking on the map always causes it to recenter. I would rather turn off click-to-recenter in general, but it's especially disconcerting when I have to bring up the menu and lose my place on the map. Is there a way to disable this behavior?


texts: Batch rectifying 200+ historic aerial images?

texts:
I have acquired historic aerial imagery for the surrounding towns, circa 1954, and the images have no georeference.

I am looking for software to help me batch-rectify the 200+ images. I have current orthophotos, DEMs, and other vector information for the area.


texts: Seeking open source alternative to File Geodatabase

texts:
I would like to use geopandas to work on a large dataset which is provided in Esri's file geodatabase format. I would like to be able to read and write this dataset, unfortunately there are still important difficulties which remain with fiona being able to read FileGDB format. Here is proof:

https://github.com/Toblerity/Fiona/issues/565
https://github.com/geopandas/geopandas/issues/868

I could of course convert everything in the FileGDB to shapefiles, but shapefiles have important limitations when compared with feature classes.
I am looking for a file format which supports the following:

Read/writeable by geopandas
Which can handle large datasets containing 10's of millions of records
Not subject to shapefile limitations
Easily convertible to and from an Esri FileGDB

Anyone know of an alternative?


texts:
I have rendered the world map using world-110m topojson. I directly used the code provided by Jakob on D3 mapping tutorial - http://www.digital-geography.com/d3-js-mapping-tutorial-1-set-initial-webmap/#.V4iGN9J97IV

<!DOCTYPE html>
<meta charset="utf-8">

<!-- Set a style for our worldshape-data -->
  <style>
  path {
    stroke: red;
    stroke-width: 0.5px;
    fill: grey;
  }
  </style>
<body>

  <!-- implementation of the hosted D3- and TopoJson js-libraries -->
  <script src="http://d3js.org/d3.v3.min.js"></script>
  <script src="http://d3js.org/topojson.v0.min.js"></script>

  <!-- map creation -->
  <script>
  // canvas resolution
  var width = 1000,
      height = 600;

  // projection-settings for mercator
  var projection = d3.geo.mercator()
      // where to center the map in degrees
      .center([0, 50 ])
      // zoomlevel
      .scale(100)
      // map-rotation
      .rotate([0,0]);

  // defines "svg" as data type and "make canvas" command
  var svg = d3.select("body").append("svg")
      .attr("width", width)
      .attr("height", height);

  // defines "path" as return of geographic features
  var path = d3.geo.path()
      .projection(projection);

  // group the svg layers
  var g = svg.append("g");

  // load data and display the map on the canvas with country geometries
  d3.json("world-110m.json", function(error, topology) {
      g.selectAll("path")
        .data(topojson.object(topology, topology.objects.countries)
            .geometries)
      .enter()
        .append("path")
        .attr("d", path)
  });

  // zoom and pan functionality
  /*var zoom = d3.behavior.zoom()
      .on("zoom",function() {
          g.attr("transform","translate("+
              d3.event.translate.join(",")+")scale("+d3.event.scale+")");
          g.selectAll("path")
              .attr("d", path.projection(projection));
    });

  svg.call(zoom)*/

  </script>
</body>
</html>

Now i want to reuse this code to render the map of India. But when I link it to India topoJSON file, only a blank SVG is created. The js console gives an error - topojson.v0.min.js:1 Uncaught TypeError: Cannot read property 'type' of undefined

I have placed world-110m.json and india.json on dropbox -
https://www.dropbox.com/sh/wrxyngyq4jdie9c/AACG2-dTzC79rRvlxlOC5poBa?dl=0


texts: Trouble in adding a print function to GeoExt MapPanel

texts:
I am trying to add print function to my GeoExt 2 web application, like:

http://geoext.github.com/geoext2/examples/printpreview/print-preview.html

first I am testing this using the same example urls, so in app.html I add:

<script type="text/javascript" src="http://demo.opengeo.org/geoserver/pdf/info.json?var=printCapabilities"></script>

then in my js file :

printProvider : Ext.create('GeoExt.data.PrintProvider', {
            method: "GET",
            capabilities: printCapabilities,
            customParams: {
                mapTitle: "GeoExt Printing Demo",
                comment: "This demo shows how to use GeoExt.PrintMapPanel"
            }
        }),

myP = Ext.create('GeoExt.panel.Map', {
            renderTo: "content",
            width: 500,
            height: 350,
            map: {
                maxExtent: new OpenLayers.Bounds(
                    143.835, -43.648,
                    148.479, -39.574
                ),
                maxResolution: 0.018140625,
                projection: "EPSG:4326",
                units: 'degrees'
            },
            layers: [new OpenLayers.Layer.WMS("Tasmania State Boundaries",
                "http://demo.opengeo.org/geoserver/wms",
                {layers: "topp:tasmania_state_boundaries"},
                {singleTile: true
                 , numZoomLevels: 8
                    })],
            center: [146.56, -41.56],
            zoom: 0,
            bbar: [{
                text: "Print...",
                handler: function(){
                    // A window with the PrintMapPanel, which we can use to adjust
                    // the print extent before creating the pdf.
                    printDialog = Ext.create('Ext.Window', {
                        title: "Print Preview",
                        layout: "fit",
                        width: 350,
                        autoHeight: true,
                        items: [{
                            xtype: "gx_printmappanel",
                            sourceMap: myP,
                            printProvider: printProvider
                                }
                            ],
                        bbar: [{
                            text: "Create PDF",
                            handler: function(){ printDialog.items.get(0).print(); }
                        }]
                    });
                    printDialog.show();
                }
            }]
        });

but I keep getting Cannot read property 'capabilities' of undefined error. I am not sure what I miss here :(


texts: Is it possible to do a query to multiple tables?

texts:
I have a point a line and a polygon layer. They have the same attribute table. Can I do one query that refers to all the attribute tables? So I wouldn't have to do every query three times. Thanks.

I'm working ith qgis, v. 2.0.1.

Edit: I joined the layers but how am I doing a query that covers all the attributes? For example I'd like to ask: Which of the attributes have the name "tree"? As I have some trees along a river, I digitized them as a line, but as I have also single trees, I also have "tree"-attributes as points. How can I do just one query for this?
Any suggestions?


texts: QGIS 3 crashes when trying to download semi-automatic classification plugin

texts:
Installed QGIS 3 in prescribed manner. (i.e. FIRST python, then GDAL, then QGIS 3).
Then installed matplotlib manually by typing "sudo pip3 install matplotlib" into terminal.
Then, (as per instructions on https://fromgistors.blogspot.com/p/plugin-installation2.html) I went to  https://www.kyngchaos.com/software/python/  to download:

NumPy 1.9.2-1
SciPy 0.16.0-1

Did not download matplotlib 1.4.3-1 because I believe I installed that when I typed "sudo pip3 install matplotlib" into terminal.

Rebooted macbook (OS 10.14.2)
Started QGIS 3 alright.
Went into manage plugins and clicked : "re-install semi-automatic classification plugin". QGIS 3 just crashed.

What's going on ?


texts: Import Shapefile to GDB. OGR2OGR FGDB Driver bug(?) on CentOS 7

texts:
I used ogr2ogr to convert shapefiles to fgdb feature classes.
I have like 10 of shapefiles so I imported the folder (input as a folder directory) and converted to FGDB
This is the code I used

ogr2ogr -progress -overwrite -f "FileGDB" $HOME/MyData/Output.gdb "$HOME/MyData/MyShapefiles/";

ogrinfo $HOME/MyData/Output.gdb seems like I have all the layers that I want to convert

However, check out the screenshot below.
This is the original shapefile

This is the one in Geodatabase

something has changed and I overlayed the two layers

Is this a bug? or is this centOS's issue??


texts: Elevation adjustment DSM using Ground control points

texts: is there any possibility to improve the vertical accuracy of a given satelite DEM/DSM by using exact measured ground control points? And if so, which software is recommended?


texts: Bounding box argument ignored when displaying raster using mapserver

texts:
I am trying to display a raster image with Mapserver using different bbox arguments in a mapserv cgi request. I have been trying to specify different bbox arguments to "crop" the raster display, but mapserver returns the full map ignoring the bbox argument.

How can I corectly specify the bounding box to display the raster?

The following example call renders the map at its full extent, although only a small part is requested

http://urbmet.webfactional.com/mapserv.cgi?MODE=map&MAP=gstack_example.map&LAYERS=HYP_50M_SR_W&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A4326&BBOX=0,0,10,10&WIDTH=128

Here is the mapfile I am using:

MAP
SHAPEPATH     "shapefiles"
IMAGETYPE     PNG
EXTENT        -180 -90 180 90
SIZE          400 200
IMAGECOLOR    255 255 255
UNITS         DD
STATUS        ON
DEBUG         ON

PROJECTION
    "init=epsg:4326"
END

LAYER
    NAME      "HYP_50M_SR_W"
    DATA      HYP_50M_SR_W/HYP_50M_SR_W.tif
    STATUS    ON
    TYPE      RASTER
END
END


texts: Skip errors during QGIS batch processing

texts: Is there an option to skip any files which result in an error when batch processing in QGIS? I'm attempting to batch process about 100 shapefiles, some of which have null geometries. Processing stops as soon as processing reaches the first bad file. Can I set it to skip these files? Ideally I'd like to be able to see which files were skipped, but this isn't essential.


texts: OpenLayers 3 for consuming an INSPIRE WMS

texts:
I am just getting into OpenLayers3 (OL3) from a GIS background - I don't have much of a programming background.

I am trying to use INSPIRE WMS layers in an OL3 app (http://www.geostore.com/environment-agency/WebStore?xml=environment-agency/xml/inspireDataWms.xml)

I have used the tutorial below as a start and I can get that working but I can't get any INSPIRE WMS layers to work.

http://astuntechnology.github.io/osgis-ol3-leaflet/ol3/README.html

I am using QGIS for the GIS work and I have checked the WMS works in that fine.

Code for the layer I have tried:

var lyr_NationalParksEngland = new ol.layer.Tile({
                    source: new ol.source.TileWMS(({
                      url: http://www.geostore.com/OGC/OGCInterface?UID%3DUDATAGOV2011%26PASSWORD%3Ddatagov2011%26INTERFACE%3DENVIRONMENT%26LC%3D400000000000000000,
                      params: {"LAYERS": "eainspire2011-wms-nationalparks_eng_inspire", "TILED": "true"},
                    })),
                    title: "National Parks (England)"
                  });


texts: How can I import 3D model in QGIS from this City Building data

texts:
I want to get the 3D model from Geoportal.  I am using the website "https://geoportal-hamburg.de/geo-online/#"
.See fig-1.  I searched for the Building and got a feature named ' Building LoD,. There I got a file called
"http://daten-hamburg.de/geographie_geologie_geobasisdaten/3d_stadtmodell_lod2/LoD2-DE_HH_2021-04-01.zip (Fig-2)"
I downloaded that data. Within the folder, there is a lot of files with Extention .gfs.
While I was trying to open it in QGIS I got an error message.  (Fig 3)

Is this possible to open 3D building data in QGIS? What's the file format.
2.How can I import 3D file. Which file format should I look for?


texts: calculation with raster and numeric values in Loop with R

texts:
I want to do some calculations with a list of rasters of lakes. Therefore I wrote a loop, which has so far worked fine when I only used rasters. Now I have calculations with  rasters and numeric vavlues and somehow this causes an error in the overlay function.
My code is:

list files

NDVI <- list.files(path =
"C:/Users/Felix/Desktop/Bachelorarbeit/Daten/Datenverarbeitung_R/L8_NDVI",
               pattern = 'NDVI.tif$', full.names=T)

define function

fun_PV <- function(x,y,z){((x-y)/(z+y))^2}

Loop

for(i in seq_along(NDVI))
{

  PV <- overlay(x = raster(NDVI[i]),
                y = minValue(raster(NDVI[i])[[i]]),
                z = maxValue(raster(NDVI[i])[[i]]),
                fun = fun_PV)

  sat = substr(NDVI[i], 73, 77)
  path = substr(NDVI[i], 78, 84)
  date = substr(NDVI[i], 85, 92)

  setwd
  ("C:/Users/Felix/Desktop/Bachelorarbeit/Daten/Datenverarbeitung_R/L8_PV")
  writeRaster(PV,filename = paste0(sat, path, date, '_PV', '.tif'))

  removeTmpFiles(h=0.1)
}

The error I get is:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘overlay’ for signature ‘"RasterLayer", "numeric"’


texts: How to convert multipoint layer to polyline layer in QGis

texts: How can I convert a multipoint layer to polylines layer in QGis. I am looking for a tool equal as the GeoWizards (convert > Multipoint to polyline) for ArcGis.


texts: Extract and sort x coordinates in groups from a column to a new one in geopandas?

texts:
I am working with geopandas

There are these columns:geo,cat and I have to make the new column called
rank.

rank should contain according to the groups with the similar cat number an increasing number according to the (about to be sorted) x coordinates which will be extracted from the geo column to a new column that will only contain the x coordinates.

  geo
  POINT (270504.6944782521 4277096.25338444)
  POINT (270838.1885699595 4278145.324786565)
  POINT (270606.3947049045 4277995.041739198)
  POINT (271508.653647932 4278548.893014569)

It will look like that

    cat   rank sorted_x
0   100   1    #here will be the westernmost point of the values where cat is 100 and will proceed to the easternmost. Then the same for the 101 values in cat and so on.
1   100   2
2   100   3
3   101   1
4   101   2
5   102   1
6   102   2
7   103   1
8   103   2
9   103   3
10  103   4
11  104   1
12  104   2

To sum up: extract x coords from the geo to another new column then group by their cat column and do the ranking on a new rank column with a pattern (for example from west to east) based on the x coordinates(which will then be sorted in order to have this output) of each group.

Here is what I have done already to achieve this result:

df['rank'] = df.groupby('cat').cumcount()+1.astype(str).str.zfill(2)

   cat   rank
0   100   01
1   100   02
2   100   03
3   101   01
4   101   02
5   102   01
6   102   02
7   103   01
8   103   02
9   103   03
10  103   04
11  104   01
12  104   02

The problem is that it doesn't do it based on the x coordinates but I think it should help to pass the logic of what I want to do.


texts: Snapping to currently digitized feature using QGIS

texts:
In QGIS 3.14 there's a new feature "Snapping to the currently digitized feature". This feature is demonstrated within a gif on the QGIS website (also shown below).
How can this feature be fully utilized?
Such as...

How can the 90-degree angles be displayed?
How can the curser position and angle be displayed?
How to properly use the "Alt + A" function?

I've been searching for a tutorial on how to use this feature, and I've had no luck.


texts:
Am using Modis Reprojection Tool (ModisTool) on MODIS LST data (MOD11A1) to
convert to GeoTiff format.

When I run it I get this error message:

Error: OpenGeoTIFFFILE: Opening Output Image File
Unable to open C:\(directory)\ ..LST_Day_1km.tif: Cannot open.

What is the correct syntax to use? Is there some way to debug this error?


texts: Python library or algorithm to generate arc geometry from three coordinate pairs?

texts: I'm trying to write code to generate fixtures for a django/geodjango project. I need to generate geometry for arcs (line string) given three lat/long pairs; start of arc, end of arc and centre of arc. I read through the geodjango documentation but this functionality doesn't seem to be available. Does anyone know if there is a Python library that provides this functionality? Or does anyone have an algorithm I can port to python? Thanks


texts: Any low-cost Virtual Private Servers able to run GeoServer?

texts:
I'm looking at virtual private server hosting options for GeoServer. At a recent conference, the most often cited  cost of such servers was around $50/month. I've found a few places that advertize VPSs for much less, some as low as $17/month. The specs appear to be adequate for installing and running GeoServer, at least for my initial experimental needs. The reason I'm looking at the lower cost options is that I'm not going to be generating any income here.

I'm curious to hear from others that may have tried this ultra-low-cost GeoServer route. How did it work out for you, and if so, which low-cost VPS providers seemed to work out the best, if any did at all?


texts: Network dataset not being built (greyed out) using ArcGIS Pro

texts:
I am trying to build a network dataset on ArcGIS Pro using a shapefile that contains polylines for different functional classification roads. I have added "length" of the polylines and "speed limit" for the roads in the shapefile attribute table. Later, I tried to build a network dataset with Time and Distance costs so, I can use it for OD matrix and routing.
But the problem is after "Building" it, the network looks greyed out and clicking on the layer does not bring Explore Network tab. After trying to create a OD Matrix based on this network dataset, it says the network is not properly built or does not show any value for TIME or Distance.
What could be the problem?


texts: Deleting lines automatically which are not connected

texts: I have a shapefile consisting of polylines. I have managed to explode these polylines (break at there vertices). There are a bunch of lines all over the shapefile which are not connected (shown in the red boxes). Is there a way to delete them automatically. What I mean is not selecting them individually because its difficult to manually locate each and every line which isn't connected.


texts: Using Esri Web AppBuilder Directions Widgets

texts:
I have a road network which I published on ArcGIS portal. The road network is working fine on ArcMap and when I published it is working fine using ArcGIS REST API directly (.../RoadNetwork/NAServer/Route/).

I am having issues while configuring the directions widget in Web AppBuilder (2.13 or before).

Travel Mode URL is giving me suggestion to add a GP Service why should I use a GP Service when I have the option to retrieve travel mode using ArcGIS rest API (..../RoadNetwork/NAServer/Route/retrieveTravelModes/) does this mean I have to publish a GP Service ?

I am not giving anything in Travel Mode URL option and When I configure the widget and run the directions I am getting following error in console.

What I am doing wrong ?


texts: Atlas Generation attribute table feature filter

texts:
I have generated a map using atlas.

I have an attribute table on the map, in this table is a column showing features. There are no records in this field, how can I turn it off / hide it so that it is not shown on the map?


texts: Is there any "proper" way of setting the viewport for QGis Time Manager?

texts:
I'm working on video visualisation for some geo data. Everything works fine for me except the viewport that the time manager plugin is using for image frames. It makes image frames size of actual map window size which is something like 1248x636 in my case now. When I resize the QGIS window the time manager plugin starts exporting image frames of different size.

Is there any "proper" way of setting the viewport so I don't have to play around with application window size? Or should I perhaps bulk edit (crop?) all the frame images before merging them with mencoder?

I'm using Quantum GIS 1.8.0 with Time Manager 0.7


texts:
I'm using ol.View option constrainResolution: true so that mouse scroll zooming has the same zooming steps as clicking in + or - zoom control. There are no custom view resolutions set, just plain XYZ map and plain view:

var openTopoLayer = new ol.layer.Tile({
  source: new ol.source.XYZ({
    url: 'https://{a-c}.tile.opentopomap.org/{z}/{x}/{y}.png',
    attributions: '&copy OpenTopoMap'
  }),
  title: 'OpenTopo',
  type: 'base'
});

var view = new ol.View({
  projection: currProjCode,
  center: ol.extent.getCenter(currProjExtent),
  minZoom: minZoom,
  maxZoom: maxZoom,
  zoom: initialZoom,
  smoothResolutionConstraint: false,
  constrainResolution: true
});

var map = new ol.Map({
  target: 'map',
  layers: [openTopoLayer]
});

Mouse wheel zooms works OK as far as resolutions are concerned, but it is not smooth. It is performed in two steps: first to half of the zoom step, then it stops and then the rest. Clicking on  + or - zoom control zooms smoothly all the way.

I tried to find related option setting in controls and interactions, but without success. The closest seemd to be maxDelta option of ol.interaction.MouseWheelZoom, but it didn't help.

EDIT: After some more testing I found out that this problem appears in version 6.2.1, but not in version 6.1.1, so it seems to be a bug.


texts:
I'm pretty comfortable with postgis, but I realized that for part of my application, having a database is pretty heavy.  I'd like to be able to take a stream of Point objects and check against an in memory collection of Polygon objects for intersection.  With postgis, this would be accomplished with ST_Contains().

The biggest reason for moving away from a database is I need real time analysis of the Points as they come in.  Right now, I collect a series of Points then do batch work on them, but I find my system struggles to keep up with the input sometimes, mostly because of the I/O associated with accessing the database rather than handling the Points in real time.

I've come across GEOS, but I'm not really sure how to get started using that library.

I'm on a Mac but also have access to a Windows machine if there is a C# solution out there.

EDIT

It's been wisely pointed out that it might be more beneficial to state my problem with my current setup before trying anything too radical.

We collect 1,000 points per second of x,y coordinates.  We map these coordinates to zones (polygons) based on which zones a point intersects with (we can have overlapping zones, and thus a point can belong to more than one zone).  Once we've established which zone a point is in and which sensor reported that point, we apply our business logic to determine if any actions need to take place.

My first approach was to use the database (Postgres 9.3/postgis 2.1) to check every point as I received it.  So I would get a reading from a sensor, figure out which sensor sent me the reading, then figure out which zones that point was in, then apply my logic.  I could check a point in about 3 miliseconds and was therefore behind the 8 ball when collecting 1,000 points/second.

My next approach was to batch these zone checkings and business logic interpretations.  Now, I grab whatever new data the sensors have produced (usually 5-10 points), figure out which zone each point is in in parallel using the database, then apply my business logic to all those points, again in parallel thanks to the database.  Then I poll for new data and repeat.

This approach seems to work well, but is a hog on the system and very fragile.  I'm barely keeping up with the data.  Someone on the postgresql IRC channel pointed out I'm seeing a lot of overhead writing points to my database just to pull them out miliseconds later to do calculations on those points.  They suggested I look at real time steam processing, and that's what I'm trying to investigate now.

I'm thinking if I can read the raw sensor data and clean it before inserting it into the database (identify which zone the point is in and which sensor sent me the reading), I can then apply my business logic in the way I currently do, I should just have more wiggle room.


texts: Chart on every marker with Mapbox Gl JS API

texts:
I am trying to move from Leaflet to Mapbox GL JS API. In my current application, I used the Leaflet Data Visualization framework (DVF) which is an awesome library. But I am afraid, is there any similar library or any solution for Mapbox API? I was searching for a similar solution but couldn't find any.

In the following example, you can see the different data visualization on markers using Leaflet DVF.

http://humangeo.github.io/leaflet-dvf/examples/html/markers.html

How can I get a similar solution for bar chart or Pie chart on markers using Mapbox GL JS API?


texts: Prevent duplicate layers in OpenLayers 3

texts:
After I created a map in OpenLayers 3 (with GeoServer), when I zoomed it out, it duplicated all layers as you can see in the following example:

http://openlayers.org/en/v3.2.1/examples/drag-rotate-and-zoom.html

I don't want the user to be able to drag the map left/right till infinity.

How can I set the extent (to the whole world) so I will not see the duplicated layers?


texts: How do I filter an image collection from properties of one feature collection and metadata of another image collection?

texts:
I have a feature collection (fc) with 29 features. Each feature has a "date" property, in the format "YYYY-MM-DD", the same format used to filter images by the filterDate() function.
var fc = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(-73.96, 40.781), {date: '2018-01-05'}),
  ee.Feature(ee.Geometry.Point(6.4806, 50.8012), {date: '2018-07-04'})
  // 29 more features
  ]);

For each feature there is an image. According to the geographical location and the date stored in the "date" property (e.g. for feature 1, I have a Sentinel-2 image that overlays its geometry ee.Geometry.Point(-73.96, 40.781) and that was acquired on a date close to the date stored in the "date" property date: '2018-01-05').
In this way I have created a collection of 29 images.
What am I looking for?
Now, based on this collection of Sentinel images, I would like to obtain a collection of Landsat images that were acquired on dates close to the acquisition dates of the Sentinel images and that overlap geographically.
I have tried to solve this problem in various ways and have not been able to get a handle on it.


texts: Horizontaly flipped/mirrored image using gdalwarp

texts:
I am trying to mask a raster using gdal_calc (and gdalmerge). To do so, I first bring them to the same extent and resolution using gdalwarp. However, whatever I do, my raster is flipped horizontally.

gdalwarp inout.tif output.tif -te -180 84 180 16 -ts 86400 16320
Warning 1: -ts values have miny >= maxy. This will result in a vertically flipped image.

Gdalinfo gives me the following:

Driver: GTiff/GeoTIFF
Files: input.tif
Size is 40321, 20161
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.004464285000012,90.004464284999997)
Pixel Size = (0.008928570000000,-0.008928570000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0044643,  90.0044643) (180d 0'16.07"W, 90d 0'16.07"N)
Lower Left  (-180.0044643, -90.0044355) (180d 0'16.07"W, 90d 0'15.97"S)
Upper Right ( 180.0044067,  90.0044643) (180d 0'15.86"E, 90d 0'16.07"N)
Lower Right ( 180.0044067, -90.0044355) (180d 0'15.86"E, 90d 0'15.97"S)
Center      (  -0.0000288,   0.0000144) (  0d 0' 0.10"W,  0d 0' 0.05"N)
Band 1 Block=40321x1 Type=Byte, ColorInterp=Gray

What is even stranger is, that the rasters align in QGIS. Is there any way to fix this?


texts: Reclassifying sub-property objects as properties in GEE

texts:
I am replicating the Hansen Global Forest Change tutorial from GEE to calculate annual forest loss by country in Africa. My code so far is:

//Country shapefiles
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');
var Africa = countries.filter(ee.Filter.eq('wld_rgn', 'Africa'));

//Hansen Global Forest Change image with multiple bands
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015').clip(Africa);
print("Hansen data", gfc2014);

var lossImage = gfc2014.select(['loss']);
var lossAreaImage = lossImage.multiply(ee.Image.pixelArea());
Map.addLayer(lossAreaImage);

var lossYear = gfc2014.select(['lossyear']);

var lossByYear = lossAreaImage.addBands(lossYear)
  .reduceRegions({
  reducer: ee.Reducer.sum()
  .group({
      groupField: 1}),
  collection: Africa,
  scale: 30
});

print("loss by country and year", lossByYear);

When I print the feature collection, I find that the forest loss data is stored not as a property, but as a sort of sub-property, under the label "groups". How do you reclassify these variables such that they export?

Original code here, but these are all GEE datasets, so the above will run fine.

Edit update: I also tried to extract the objects by year:

var loss0 = lossByYear.map(function(f) {
  var loss2000 = ee.Feature(ee.List(f.get('groups')).get(0))
      .get('sum');
  return f.set('loss2000', loss2000);
});

And get the error: FeatureCollection (Error)
Error in map(ID=00000000000000000000):
Element.get, argument 'object': Invalid type. Expected: Element. Actual: Dictionary.

I also tried converting the objects as follows:

var loss = ee.FeatureCollection(ee.List(lossByYear.get('groups'))
  .map(function(el) {
    var d = ee.Dictionary(el);
    return ee.Feature(null, { yearString: ee.Number(d.get('group')).format("20%02d"),
                              lossArea: d.get('sum'),
                              yearNumber: ee.Number(d.get('group')).add(2000)
    });
  }));

But get the error: FeatureCollection (Error)
List.map: Parameter 'list' is required.


texts: Repairing broken links for new paths of layers to shapefiles using ArcPy?

texts:
I have an .mxd with broken link layers and each layer has a different new path.
I know the python script for repairing the broken links of the layers, but I want to put all the new paths of the layers and then let python do and finish all the process. The script that I have does it one at a time. I know if all the shapefiles are in one path then it's done by performing it once. but what if they al have different paths?

To clarify, let's say I have 3 layers with broken links and each one has a different new path now. I'd like to have a GUI with 3 inputs for the new paths altogether or let script do the process one by one and each time asking for the next path to input.


texts: QGIS How to multiply (repeat n times) a point by field

texts:
I have a point layer in QGIS that has a field named "Amount" that stores an integer value for every point in the layer.

I want to create N points in the exact place of the feature that has N as value in the amount "field". I mean, if a point has in the amount field 16, I need to create 15 more points in that exact lat long (1+15 so that there are 16 in total) of that point.

After this process, the sum of "Amount" values in the starting layer should be equal to the sum of points in the new layer.

I need it to create a heatmap with the new layer. What is the easyest way of achiving this?

Regards


texts: Date field formatting as number in ESRI popup window

texts:
We are using the Park and Recreation Locator from ArcGIS for State Government and trying to put a date field into the popup info window.  The field is pulling the correct data, but it's formatting it as a number (JavaScript milliseconds) instead of a date.  There is a variable in the config.js file called FormatDateAs which has been set as a date according to ESRI's site "MMM dd, yyyy", but has no discernible effect.  From experience with some of ESRI's other templates, I tried adding in fields to the InfoPopupFieldsCollection such as isDate or variations on Format but none worked.  ESRI's site has no additional documentation for formatting a date in this template.

We need the field to remain a date type in the database for future changes we have planned, so simply changing the fields to strings containing the pre-formatted dates won't work.  Does anyone know how we can get this field to format correctly in the code?


texts: Calculating Field Chosen by arcpy.GetParameterAsText(1) gives ERROR 000539?

texts:
I add a parameter CampoConfere = arcpy.GetParameterAsText(1) which get a field from Tabela_Base to calculate in the next steps from the code. The Scripts works well until the Calculate Field Step when the prompt says that "CampoConfere" is not valid. Why Python do not calculate this parameter?

import arcpy

Tabela_Base = arcpy.GetParameterAsText(0)
CampoConfere = arcpy.GetParameterAsText(1)
arcpy.AddMessage('Tabela de Entrada: '+Tabela_Base)
arcpy.AddMessage('Campo: '+CampoConfere)
fields = arcpy.ListFields(Tabela_Base)
ValidaCampo = sum(field.name == CampoConfere and field.type == "Single" for field in fields)

if not ValidaCampo:

    Nome_Campo = "CHAVEIBGE2"
    Tipo_Campo = "FLOAT"
    Campo_Calculo = "CHAVEIBGE2"
    Saida1 = Tabela_Base

arcpy.AddField_management(Tabela_Base, Nome_Campo, Tipo_Campo, "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.CalculateField_management(Saida1, Campo_Calculo, "!CampoConfere!", "PYTHON_9.3", "")

Them I get the message error:

  Executing: Script F:\AgroDB.gdb\IBGE\Brasil_Municipios COD_IBGE
  Start Time: Fri Jun 23 08:46:05 2017
  Running script Script...
  Tabela de Entrada: F:\AgroDB.gdb\IBGE\Brasil_Municipios
  Tabela de Entrada: COD_IBGE
  Failed script Script...

  Traceback (most recent call last):
   File "F:\AgroDB.gdb\Modelo1.py", line 18, in
     arcpy.CalculateField_management(Saida1, Campo_Calculo, "!CampoConfere!", >"PYTHON_9.3", "")
   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", >line 3453, in CalculateField
     raise e
  ExecuteError: ERROR 000539: Invalid field CampoConfere
  Failed to execute (CalculateField).

  Failed to execute (Script).
  Failed at Fri Jun 23 08:46:06 2017 (Elapsed Time: 0,74 seconds)


texts: Updating table with GeometryFromText of PostGIS

texts:
I'm trying to update a column public_geometry in a table  by creating a point from longitude and latitude cols from the same table. The database is PostgreSQL with Pgadmin as program to handle the databasse:

 UPDATE sirene_test.etablisssement_actif
SET geom=ST_GeometryFromText('POINT('|| latitude  ||' '|| longitude ||')', 4326);

I got this Error :

ERREUR:  parse error - invalid geometry
HINT:  "POINT( )" <-- parse error at position 9 within geometry
********** Erreur **********

ERREUR: parse error - invalid geometry
État SQL :XX000
Astuce : "POINT( )" <-- parse error at position 9 within geometry


texts: I have an Oracle spatial layer that is built on top of a multi table view. When I make a new feature in QGIS, I get an error when I attempt to the save the edit. When I look in the database and refresh the QGIS view, the new features are present, so the error is being thrown even though the operation is successful. The errors says could not retrieve feature id -2 and has a sql statement where it looks like QGIS is quering by rowid. I was wondering if anyone elese has seen something similar to this or if it is possible to have qgis not use Oracle rowid as a check.


texts:
I am trying to extract elevation values for each point in a point .shp from a DEM. Both files align and are in the same CRS - 4326. The tool runs but I receive an error message and the output file is completely empty. Does anyone have any tips on how to get a correct result?

This error message is repeated many times:

  Feature write errors: Feature creation error (OGR error: Attempt to write non-point (MULTIPOINT) geometry to point shapefile.)


texts: Adding points to map with Leaflet using .forEach function

texts:
Trying to add an array of points to a map with a Leaflet with .forEach.

I thought this should work, but alas, it does not:

var jsonData = [
  {
    "LNG":-75.16472099534634,
    "LAT":40.034508187448395,
    "NAME": "Resto1"
  },
  {
    "LNG":-75.17898120530263,
    "LAT":40.06738292573648,
    "NAME": "Resto2"
  }]

var addMarkers = function(point) {
      L.marker([point.LAT, point.LNG])
        .addTo(map)
        .bindPopup(point.NAME);
};

jsonData.forEach(addMarkers);

I understand I can do it like this:

var addMarkers = function(map) {
  for (var point in jsonData) {
     L.marker([jsonData[point].LAT, jsonData[point].LNG])
        .addTo(map)

But was hoping to do it more functionally. Any hints?


texts: Geoserver not showing compressed rasters?

texts:
beginner here.  I'm practicing loading a raster file into Geoserver, and loading/viewing it works fine with my uncompressed (1.2GB) tif (seems I can use either GeoTiff or Worldfile options just fine), but it will simply not display compressed tifs using the same workflow.

A bit more detail of what I've done:

I merged some landsat images into a large GeoTiff via gdalwarp (1.2GB file)
I use listgeo command to create a tfw file associated with the tif
I edited the tif in GIMP to adjust some color settings, etc, which loses the georeference data
After GIMP edits complete, I use gdal_edit.py to reassign srs (EPSG:3857) and ensure tfw has same name as tif file.
At this point, I can then load the large tif/tfw files (1.2GB) onto my Geoserver as a store and then publish as layer, and I can see/zoom,etc the raster just fine in layer preview (openlayers).
However, if I use any compression ("gdal_translate -co COMPRESS=LZW -co TILED=YES in.tif out.tif" or COMPRESS=JPEG) and upload this smaller (~200MB) file (and associated tfw) to geoserver as a store and publish as layer (same procedure as uncompressed tif), it simply does not show when I try to view in layer preview (open layers) as the larger one does.  The window's coordinate area seems fine, but it is just a blank box with no visible raster!

Does anybody have any ideas why I might be experiencing this?

I figured out I was missing a step in the re-georeferencing (I used gdal_translate's -a_ullr to add in upper right, lower left, etc coords back into the tif).  So I can now get the properly georeferenced compressed tif to show, so I'm assuming it had something to do with that.  I still can't get the non-properly georeferenced compressed version to show on Geoserver, though the non-properly georeferenced uncompressed one does just fine still.  This is GeoServer v2.9.0. I'll consider this matter resolved for now.


texts: Merging contour lines along tile boundaries

texts: I am working with the CANVEC data set. The data in canvec is broken up in to mapsheets. eg. 116B04. I have data covering the entire province of the Yukon. I have merged all the data in to one shapefile but I would like to take this further by creating one feature for each line of elevation within the entire dataset. I would like to do this so that I can display things a little quicker when I load it in to QGIS but dont want to turn it in to a raster so that I can switch labelling on and off to view elevation labels. I am thinking that creating one layer out of the entire dataset and one attribute entry for each elevation. eg. 1500ft, 2500ft, 3000ft. ect. I am using UBUNTU 13.10, QGIS, GRASS, POSTGIS. I would also like to do something similar for polygons that are split into parts because of the map tiles.


texts:
I am very new to Python and my first attempt hasn't gone too well.  I'm trying to overwrite a field that is a spatial join of a couple of features in order to, in the end, generate a beginning and ending street intersections listing.
Here is a link to my table.  The field I'm trying to overwrite is the IntersectionList field.

As you can see, the IntersectionList has the SECONDARYN field names included due to the spatial join. All I need now, are those streets that are not in the SECONDARYN field in order to create a beginning/ending street listing.
Here is a visual example of what I'm trying to accomplish.

This is the code I tried to use via the Calculate Field tool:
fixUpIntersections(!IntersectionList!, !FULLNAME!)

def fixUpIntersections(listAsString, streetname):
names = listAsString.split(" & ")
uniqueNames = list(set(names))
if len(uniqueNames) == 1:
    return uniqueNames[0]
    uniqueNames.remove(streetname)
return ' & '.join(sorted(uniqueNames))

I thought with the error I received, it meant that maybe I should !FULLNAME! to !SECONDARYN!, in order to get rid of those streets in the IntersectionList but when I tried that, it just deleted all data in the IntersectionList. Running the tool with the code exactly as shows above without changing !FULLNAME!, it resulted in the following error: ERROR 000539: Invalid field FULLNAME

I've dug around the site here but haven't been able to figure out my problem. I have a long way to go to learn Python.


texts: Measuring the similarity of two polygons in PostGIS

texts:
I want to measure the similarity of two polygons. I want to use the average distance of the boundaries as a measure of the shape similarity.

How can I do this in PostGIS?


texts: Using exterior orientation parameters to calculate 3D similarity transform

texts: Given a solution from relative orientation (X,Y,Z, kappa, phi, omega) for two images (model coordinates) and exterior orientation parameters for two images from GPS/IMU (X,Y,Z,kappa,phi,omega) (object coordinates), is it possible to calculate the parameters from a 3D similarity transform relating model coordinates and object coordinates?


texts: Calculate angle of view from a point based on DEM using QGIS

texts:
I'm looking for a way to calculate the angle of view from a point on a DEM, in particular from the highest point of a terrain relief (corresponding to the red dot in this image below).

To better explain the case: I have a cell with an altitude of 237 m asl. I'd like to know the angle of "free" view from the observation point when it's not blocked from the surrounding relieves, setting a maximum distance over which the view is considered free. Whenever the DEM exceeds that altitude, the view from the observation point is blocked by the higher relieves.
Here's a graphical example of what I mean.

I've tried several times with the r.viewshed algorithm in QGIS. It let me choose a maximum distance, but it doesn't calculate the angle of view and it gives other types of outputs.
Is there a way to obtain this information?


texts: Unable to save attribute table after creating new attribute in QGIS?

texts:
I have little to no GIS experience. My workplace has no GIS capabilities.

I am using QGIS 3.0.2 (64 bit) on Windows and I created a new attribute for a vector layer (shapefile) through the field calculator by concatenating 5 other attributes. I first tried to do it with the attribute table open (which froze my machine), then tried it without opening the attribute table and it worked.

However, when I disable edit mode and save the attribute table, QGIS just keeps thinking and thinking. I left it running overnight but it was still going in the morning and I ended up shutting it down.

There are approximately 1.9 million entries in the attribute table to be saved and I understand it may take a while but should it take that long to save or am I missing something?

Or is there a better way for me to create a new attribute based on existing attributes so that I can join it to other layers in the future (join based on this new attribute)?

I am concatenating 1 string field, 4 integer fields & hyphens into a 14 character (max) string field (i.e. NW-15-25-15W5 or SE-24-112-10W4).


texts: ValueError: too many values to unpack (expected 2)

texts:
I am trying to follow this tutorial: land cover fastai
And while executing the workflow, I got this error:

This is the code where the error occurs:

def png_gen(patch, inference=False):
    """Save RGB, spectral info and labeled images for the coming deep learning

    Args:
        patch: Saved eopatches for deep learning LULC training and prediction
    Return:
        None: images in the PNG.
    """

    patch_data = EOPatch.load(patch)
    patch_dir, patch_id = patch.split("/")
    if inference == True:
        inds_path = "inds_inference"

        if not op.isdir(inds_path):
            makedirs(inds_path)
        for i in range(len(bands)):
            #get NDVI, NDWI and NDBI from the bands
            inds = color_scale(bands[i][..., [-3, -2, -1]]).astype("uint8")
            Image.fromarray(inds).save(op.join(inds_path, patch+'_'+str(i)+'.png'))

    else:
        lulc = patch_data.mask_timeless['LULC']
        lulc_path = "lulc_all"
        rgb_path = "rgb_all"
        inds_path = "inds_all"
        if not op.isdir(lulc_path):
            makedirs(lulc_path)
        if not op.isdir(rgb_path):
            makedirs(rgb_path)
        if not op.isdir(inds_path):
            makedirs(inds_path)
        for i in range(len(bands)):
            # Switch R, G, B band index from bands
            rgb = color_scale(bands[i][..., [2, 1, 0]]).astype("uint8")
            inds = color_scale(bands[i][..., [-3, -2, -1]]).astype("uint8")
            Image.fromarray(rgb).save(op.join(rgb_path, "{}_{}.png".format(patch_id, str(i))))
            Image.fromarray(inds).save(op.join(inds_path, "{}_{}.png".format(patch_id, str(i))))
            Image.fromarray(lulc).save(op.join(lulc_path, "{}_{}.png".format(patch_id, str(i))))

And this is the execution code:

#download eopatch for the desired AOI and covert the numpy array into RGB, spectral info (NDVI, NDWI, NDBI)
# and training lulc label into PNG and saved them under "rgb_all", "inds_all" and "lulc_all"
pbar = tqdm(total=len(patchIDs))
for idx, bbox in enumerate(tile_list[patchIDs]):
    # define additional parameters of the workflow
    extra_param = {
        add_data: {'bbox': bbox, 'time_interval': time_interval},
        save_s2: {'eopatch_folder': 'eopatch_{}'.format(idx)}
    }

    workflow.execute(extra_param)
    print("eopatch {} has been processed!".format(idx))
    # png_gen('eopatch_{}'.format(idx))
    png_gen(op.join(path_out_sampled, 'eopatch_{}'.format(idx)), inference=False)
    print("Saving RGB, LULC and Inds PNG to eopatch_{} for the coming ML pipeline".format(idx))
    # del 'eopatch_{}'.format(idx)
    shutil.rmtree(op.join(path_out_sampled, 'eopatch_{}'.format(idx)), ignore_errors=True)
    print("eopatch_{} deleted!".format(idx))
    pbar.update(1)

Instead of  '/' I tried  '_' in this line  patch_dir, patch_id = patch.split("/") because the downloaded patches have this format "eopatch_idx", but I got the same error.

This where the downloaded data are stored:

path_out_sampled = './eopatches_sampled_Slovenia/'
if not op.isdir(path_out_sampled):
    makedirs(path_out_sampled)
save_s2 = SaveToDisk(path_out_sampled, overwrite_permission=OverwritePermission.OVERWRITE_PATCH)


texts: Measuring distance between two coordinate pairs taken from NMEA file?

texts:
In order to measure momentary speed, first I need to measure the distance between two pairs of coordinates. That's where my lack of knowledge stops me.

I have example data in NMEA format:

$GPGGA,002246,3918.5986,N,11954.1140,W,1,08,2.3,2697.9,M,-22.4,M,,*4C
(...)
$GPGGA,002248,3918.5887,N,11954.1395,W,1,08,2.4,2696.5,M,-22.4,M,,*42

I figured out that data above gives me following information, needed to measure the distance:

Coordinates: 39 deg 18.5986' N, 119 deg 54.1140' W
Altitude (above sea level): 2697.9 m
Height of geoid above WGS84 elipsoid: -22.4 m

Coordinates: 39 deg 18.5887' N, 119 deg 54.1395' W
Altitude (above sea level): 2696.5 m
Height of geoid above WGS84 elipsoid: -22.4 m

Angular distance between this two points is: latitude -0,0099', longitude 0.0255'.
The difference between altitudes is -1,4 m.

And there I'm stuck. And my question is - how to measure the distance between those two points in (kilo)meters?

I think it's worth repeating that the distance would be very small, because it's needed to calculate a momentary speed. So I think the data emphasized  above (the height of a geoide) could improve precision.


texts: How do I run SAGA algorithms in PyCharm?

texts:
Is there anyway I can input a raster and number and then use the SAGA algorithm "Fill Sinks (Wang and Liu)" on the raster through an IDE such as PyCharm?

How would I go about calling the algorithm?

What I want to do is write code in Python that inputs a raster and then outputs it by calling a function from the SAGA toolkit.


texts:
I have 1 vector layer (counties) and 1 raster layer (population density). I have tried merging these two, so that I can examine population density within a county, but I have had no luck. It just gives me columns with NULL.

I did the standard thing, which is to calculate Zonal Statistics. My raster layer is an asc file, and I have tried converting it to a tif. Nothing has worked.

I got my population data from the History Database of the Global Environment (HYDE).

Okay, it could be because my CRS projections aren't the same for the vector and raster layers. When I set both the same, however, the map completely screws out of whack, with the vector layer being off-map completely! What do I do?


texts: Using pysal to map location of coordinates?

texts:
I am trying to map which county a set of coordinates I have are in. I want the code to go through each coordinate, check which county it is in, and then print the corresponding county. However, currently when ran the cursor in idle just flashes and nothing appears to happen. I'm not sure if the code is correctly checking the counties or whether it can get the the names of the counties from the shapefile?

import pysal as ps
from ast import literal_eval

COUNTIES = '/Users/JoshuaHawley/Desktop/counties/English Ceremonial Counties.shp'
shps = ps.open(COUNTIES)

counties = shps.read()

f = open('/Users/JoshuaHawley/sundayCoordinates.txt').read()
seq = literal_eval(f.replace("\n", ","))

coordinates = seq[0:]
for line in coordinates:
    points = coordinates[:2]

NotFound = True

while NotFound == True:
    if coordinates in counties:
        NotFound = False
        name_of_county = counties.getNamePolygon()
        print(name_of_county)


texts: Carrying attribution from several points to single underlaying polygon to create copies of polygon with unique attribution using ArcGIS Pro

texts:
I am using ArcGIS Pro to manipulate some property parcels data. I have a parcels polygon layer with the parcel geometry, and a points (PINs) layer with the property attribution information. Most parcels have a single PIN, so I can easily join the attribution from the PIN to the parcel using a one-to-one spatial join.
For those parcels that have several PINs within them, however, I'd like to create one copy of the underlaying parcel for each PIN it contains. I'd like each copy of the parcel to be given the attribution of one of the PINs. I've tried running a spatial join using the one-to-many join operation, but this just creates duplicate copies of the polygon with attribution from the first PIN only. It seems like this should be fairly easy but I haven't been able to find a good solution to this problem.


texts: How to scrape polygon data from web map?

texts:
for each value of parcel no: in http://www.ma-investment.gov.bh/bps_ppd_en/ the map generates polygons.

I want to extract xy values corresponding to polygon shapes which i can convert to longitude,latitude using gpoint.php


texts: Display of image in GEE (Sentinel-2) does not display the masked values using SCL band


texts: Creating dictionary of link and its adjacent links (having common vertex) in ArcPy

texts:
How do I create a dictionary that has the linkID of a link as key and linkID of touching links as values?

I have used this code in ArcPy but it takes a very long time.

arcpy.CopyFeatures_management(link,link1)
arcpy.MakeFeatureLayer_management(link,"link")
arcpy.MakeFeatureLayer_management(link,"link1")

fid_fid_dict = {} #the dictionary that has linkID and its connecting linkIDs

with arcpy.da.SearchCursor(link, ["FID"]) as curs:
    for row in curs:
        arcpy.SelectLayerByAttribute_management("link", "NEW_SELECTION", '"FID" = {}'.format(row[0])) #selects a link
        arcpy.SelectLayerByLocation_management('link1', "INTERSECT", "link", "", "NEW_SELECTION") #selects touching links
        fid_fid_dict[row[0]] = [_row_[0] for _row_ in arcpy.da.SearchCursor("link1", ['FID'])] #gets the IDs of touching links

The final dictionary should be:

dict = {linkID:[linkID1, linkID2, linkID3, linkID4]}


texts: Calculating distances with SciPy and subsequently join identifiers

texts:
I have lots of start points and a handful of destination points. Using 'Point Distance' from the toolbox takes ages to calculate and write the results. Hence, I was looking for alternatives and went the Numpy, SciPy way.

Here how I calculate the distances, and it only takes some seconds if at all.:

import arcpy
import scipy.spatial
arcpy.env.workspace = 'Drive:/Path/To/Database.gdb'
start = arcpy.da.FeatureClassToNumPyArray('Foo', ['Bar', 'SHAPE@XY'])
dest = arcpy.da.FeatureClassToNumPyArray('John', ['Doe', 'SHAPE@XY'])
arr_dist = scipy.spatial.distance.cdist(start['SHAPE@XY'], dest['SHAPE@XY'], 'euclidean')

That gives me in my case a matrix with the dimension of (670911, 39). That means, each row is representing a feature from start, while each column holds the distance for the particular point from the start FeatureClass to all destination points. The order of the destination follows the order of OIDs.

My question is: How do I get the information from 'Bar' and 'Doe' to the arr_dist and get it exported into a table in my database.

The arcpy.da.NumPyArrayToTable function requires a structured array and the dtype for every column, which seems to be pretty cumbersome considering the number of columns I have.

[EDIT]

I tried to append the information from 'Doe' to arr_dist using the following line:

import numpy.lib.recfunctions
arr_doe = numpy.lib.recfunctions.append_fields(arr_dist.transpose(), 'doe', dest['Doe'])

That has the reproducible effect that ArcGIS (10.4) crashes. When running it in IDLE a 'MemoryError' is returned.


texts: Using newline/page break tag in ArcGIS Desktop Text_Element in arcpy.mapping?

texts:
I want to update a Text_Element with a newline or page break using ArcGIS Desktop 10.0.

Does anyone know the tag or syntax?


texts: "\"-character is not recognized in field calculator string (URL field)

texts:
I have this irritating problem with "\" symbol in field calculator in QGIS. I'm trying to insert an URL field in order to open feature specific pictures later as an action, but as I insert the symbol in question, it is not recognized. See below. Preview will be the result if I hit OK.

Column's data type is text and layer is UTF8 encoded. I'm using QGIS 3.2.0 on Windows 10.

It's a small issue, but I can't wrap my head around it.


texts: Save as GeoPackage Layer Options - use of Description and Identifier

texts:
When saving a layer as GeoPackage in QGIS 3, what is the purpose of the Layer Options DESCRIPTION and IDENTIFIER?

I'm looking at how we can manage our working data better than the current mess of shape files.  When you hover over DESCRIPTION the tooltip says "Human-readable description for the layer content".  I thought this might useful to as quick way to provide some documentation/metadata as to the purpose of each layer.
However the DESCRIPTION and IDENTIFIER don't seem to be exposed anywhere that I can find. When looking at the geopackage in the Browser or DB Manager or the add layer dialog I see only the table name, no DESCRIPTION or IDENTIFIER.  Am I looking in the wrong place?


texts: Leaflet event onclick zoom in and out (JavaScript)

texts:
I just want to make onclick event in zoom in and zoom out. How can I do it?
I just confused to make it because zoom in and zoom out is something that generated by library.

Should I make new function?
How I make it possible?


texts: OpenLayers 3: Change drawing style, but keep style of previously drawn features

texts:
I created a vector and a drawing interaction:

// Create feature collection and vector
var featureCollection = new ol.Collection();
var vectorSource = new ol.source.Vector({
    features: featureCollection
});
featureOverlay = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'rgba(255, 50, 0, 0.5'
        })
    })
});
featureOverlay.setMap(map);

// create drawing interaction
drawingInteraction = new ol.interaction.Draw({
    features: featureCollection,
    type: 'Polygon'
});

map.addInteraction(drawingInteraction);

Now I would like to let a user change the drawing style, which, when called separately, can be done like this:

var s = new ol.style.Style({
            fill: new ol.style.Fill({
                color: 'rgba(0, 255, 0, 0.5'
            })
});

featureOverlay.setStyle(s);

This works fine, but it changes the style of all the drawn polygons. However, I would like to change the new ones only and keep the styling of the existing ones.

For example, I draw my initial red polygons:

Then the user switches styles by firing the setStyle() method above, but then all polygons turn green, not just the one I will set after firing the method and drawing a new polygon:


texts: Postgis 2.1: ERROR: Shell is not a line - SQL state: XX000

texts:
My current setup is PostgreSQL 9.3.6 & PostGIS 2.1.5

I have a table with contour lines (multilinestrings) that I wish to group by the attribute 'elevation', and make polygons out of them. The lines, hovewer, were derived from different sheets, so it might be that one of these lines suddently ends exactly where another (with the same elevation) begins. This is an example by selecting one of those lines in QGIS:

I get a strange error message:
ERROR: Shell is not a line - SQL state: XX000
when I try to run this query:

select st_multi(st_makepolygon(st_union(the_geom)))
from import.contour
group by elevation

however if I run

select distinct geometrytype(the_geom), st_isvalid(the_geom)
from import.contour;

I get as result:

|  geometrytype  | st_isvalid |
|----------------|------------|
| MULTILINESTRING|     t      |

So it seems to be everything ok for me,...

What is wrong here?


texts: Can OpenStreetMap display train stations?

texts:
I do not have a car and I move around on a bicycle. I thus often need to know what the closest train station is to a location I need to go to.

I notice that OpenStreetMap (OSM) displays airports in an obvious way, such as here.

Is there a way to have OSM (or Google Maps) display train stations in a more obvious way than having to zoom in enough on each city in an area?


texts: Using CopyRaster in ArcGIS ModelBuilder?

texts: I need to remove white background around a pack of rasters. When I try Copy Raster from toolbox it works perfectly, but when I try batch processing or iterator in ModelBuilder, I have rasters with black background. It seems like one of parameters (ignore background value) don't used in model.


texts: Selecting features of a shapefile using attributes of another shapefile

texts:
A newbie here trying to learn QGIS for work.

How can I select specific features of a polyline with the same values as the attributes of the buffer polygons I generated?

I have a layer of polylines with numerical values as one of their attributes.

I was able to generate buffer polygons using the midpoint of each line and used the join by attributes function to determine the maximum values of the attributes of the lines within each buffer.

I now need to select the lines that contain the said maximum values for each buffer zone. Since I have a lot of lines to work with, manually finding each one is not an option.


texts: Extract .osm data appears in axonometric on QGIS

texts:
I'm working with QGIS 2.18.

I extract .osm data from the site.
name : vernon.osm
Then i transformed it in a data base file : vernon.osm.db
Finally i used spatialite to make it appearing on QGIS. polylines and polygons appear but in axonometric (like in the picture)

The main point of all this is to export the datas in .dxf to use on autocad.

In axonometric it's not exportable.


texts:
I am currently writing some code to do multilateration using TDOA and currently struggling and wondering if anyone has any reference code that can be used as a test case with some sample data that I have!  At the moment I am unsure if it is the code I have written or whether it is the geometry of the points that is off.

One example group of packets that I have received is

station1 = [50.919149572907784,-1.005043823556842,90.729];
station2 = [50.8304,-1.21956,0];
station3 = [51.0569,-1.26224,100];

stations = [station2,station1,station3];
timestamps = [28704.175275109, 28704.175334203, 28704.175334875];

My code currently locates this to 50.910447339754654,-1.2347406519401736 which I am 99% sure is wrong!  Does anyone have any code that can verify this at all?

One bit of code I am testing against and get very similar results is https://github.com/bistromath/gr-air-modes/blob/master/python/mlat.py


texts: How to display large geospatial data in a jupyter notebook with style and tooltip?

texts:
I have a large amount of data stored on a postgreSQL database (RDS) with the PostGIS extension enabled. We are currently doing quality assessment of the data and therefore I would like to inspect/share results with my colleagues overseas. I find querying the database and parsing the result to pandas and geopandas in a Jupyter Notebook (Google Colab) extremely helpful. Connecting to the database using ArcMap or QGIS is less ideal since we need a lot of charting functionality (time series).

The approach I took so far is to add popup functionality and style to a geodataframe using folium. results for a sample dataset can be found here

A caveat is that this method does not scale very well and fails on my dataset due to the number of geometries. Data and visualization happens client-side and the method fails with complex geometries are large geodataframes. I am looking into alternatives with the following requirements:

interactive (zoom,pan)
popup (tooltip) with attribute (property) values
color / style
handle 16000 geometries, large number of vertices (200MB)
ideally server-side or serverless design.
ideally open source.

So far I've looked into

cartoframes:
does exactly what I want except that tooltips are not supported yet.

Other packages that I will look into:
 1. Arcgis API for python.
 2. Mapbox API for python.
 3. mapboxgl-jupyter


texts:
I am trying to display a network in QGIS. I have network line coordinates in a .csv file, stored in 4 columns: From_X, From_Y, To_X, To_Y.

I would like to import this file into QGIS but the option Layer> Add Layer from a Delimited Text File doesn't allow me to select more than two coordinates: X field and Y field.

Is there any option that I am missing? I have read some other questions talking about WKT, do I need my data in this format?


texts: Map point class diversity on raster efficiently in R

texts:
I have a point data set of roughly 2'000 features that represent User GPS locations. I'd like to calculate the User diversity per pixel. This is a (very slow) Pseudocode:

library(raster)
library(sp)

pts <- SpatialPointsDataFrame(someQuery)
rast <- RasterLayer(someQuery) # Defines empty raster pixels

for (pixel in pixels) {
    1. Get point subset pts_sub within distance d from pixel center
    2. Get user_diversity of pts_sub (eg nUsers / nPts)
    3. Set pixel value to user_diversity
}

The following image illustrates the situation. For the upper left pixel, all pts within the blue circle have to be subset and the diversity index is calculated from these.

The procedure has to be applied to many more datasets, therefore I need a faster way to do this. The "calc"-function of the raster package allows custom functions on raster, but doesn't allow to access the pixel coordinates. Is there a fast way to apply such a function to pixels of a raster?


texts: Why the "shape tools" won't work with a different CRS?

texts:
I just used a tool from the "Shape Tools" plugin but the results were nonsense. The CRS of the point layer and my project CRS were  EPSG:4326 - WGS 84.

Here is an image of the result: The correct result should be an ellipse around the point.

I realized if I change both the project CRS and point layer CRS to EPSG:25832 - ETRS92/UTMzone 32N. The tools of the plugin work perfectly fine.

That's the correct ellipse.

Any explanaition why this is the case?


texts:
I was able to login as

http://<domain_name>:3000

Done some import jobs, publish the map, copied the API link to a sample HTML page. Able to view the whole base map and CartoDB layer from local. But I tried from other network, the CartoDB layer seems to be not seen. Checked the console, but no errors where logged. Edited the development.js from the Windshaft directory and changed the localhost instances to domain name, and also the app_config.yml files.


texts:
I must process an  Excel file data  but  the format of my longitude and latitude coordinate are: Lat: -22 to -34 and Long: 9 to 18 in decimal degree. So i want to use Add XY tool into Arcgis10 to put my coordinate data, but i have meet some difficulties because two latitude and two longitudes.
How to put the variation of latitude and longitude coordinates with Arcgis tool?


texts: How to add multiple relative hyperlinks to single point in ArcGIS for Desktop with Basic level license?

texts:
It may be a very obvious one, but how do I add more than one relative hyperlink to a point feature class in ArcGIS for Desktop with a Basic level license?

The Basic level license eliminates using Attachments being a possibility.


texts:
I am trying to create a spatial index for some shapefiles using the Qgis API. I have followed the steps described in Nathan Woodrow's blog (https://nathanw.net/2013/01/04/using-a-qgis-spatial-index-to-speed-up-your-code/) and it works great:

layer = QgsVectorLayer(path, name, 'ogr')
idx = QgsSpatialIndex()
all_features = layer.getFeatures()
map(idx.insertFeature, all_features)

My problem is that I don't have any resulting file (.qix ? .sbn ? .sbx ?)

How can I save this spatial index so that future users of the shapefiles can benefit from it ?


texts: Postgres/PostGIS: Dissolve multiple fields with COUNT (Point shape)

texts:
I have loaded point .shp feature classes into Postgres 9.3 db, and enabled postgis,plpgsql,postgis_topology extensions.

I'm looking to dissolve multiple fields within a single table 'TABLE1' and add a COUNT or similar field to the tabular output.  The feature class has more than 12Million points, so Esri understandably chokes, and QGIS only allows me to dissolve based on one field at a time.

Start table:

The summary table I am trying to generate:

So, I'm looking to summarize across all 3 fields (ORG, PROTOCOL, MAP) and then also generate the 'Cnt_Map' field (or similar) which denotes the COUNT of each MAP.

Is there a specific SQL command or other query I can use to run this type of report?

Grazie


texts:
What is the delta between two lat/lon pairs called? Delta degrees? Arcs?

Is it some quantity of distance but what is its unit of measurement?

And related: Can we use pythagoras to calculate this distance?

eg. d = (0, 3) and (4, 0) = 5?


texts: QGIS dissolve/merge features by all attributes

texts: I have a vector layer with thousands of features and over ten columns (attributes). I would like to dissolve/merge features with identical attribute combinations and keep the attributes. I think it could be "singleparts to multiparts" but I can´t find this in the toolbox. In dissolve by attribute it is only possible for three attributes. I remember that I did it very easy one day...


texts: Intersecting polygon of certain elevation with TIN surface?

texts:
I want to generate a flood plain for a river in ArcMap.  Basically, this means that I want to create polygons with certain elevations above the drainage point elevation (870 m).

For instance, I would like to make a polygon that has an elevation of 872 m, intersect that polygon with the surface tin which would give me an idea of the flood plain if the water level was two meters above the drainage point.


texts: Exporting records from an attribute table and placing them in a txt file

texts:
I am trying desperately to grab records in a shapefile's attribute table and copy them, with additional formatting to a txt file.  Below is the code I have developed thus far:

I have used the following post (Using SearchCursor to access & export values in a raster table) for guidance but for some reason python has issues with my declared output file, f.

# Import arcpy and the SeachCursor function

import arcpy
# Find the target shapefile in memory
shapefile = "C:/Users/user/Dropbox/CE 691/GI/GISData/ExampleEPANETNetwork.mdb/Nodes"

# Create the Text File to write the Indicated Values
f = open("C:/Users/user/Desktop/test.txt", "w")

# Use SearchCursor to Identify the necessary fields
find_junction = arcpy.SearchCursor(shapefile, "", "",  "Elevation;Demand;X_Coord;Y_Coord","")

# Convert Found Field to Text and Place in Text File
for row in find_junction:
    get_elevation = row.getValue("Elevation")
    get_demand = row.getValue("Demand")
    get_x_coord = row.getValue("X_Coord")
    get_y_coord = row.getValue("Y_Coord")
    print get_elevation, get_demand, get_x_coord, get_y_coord
    f.write(str(get_elevation + "    " + get_demand + "    " + get_x_coord + "    " + get_y_coord)

# Close the output text file
f.close()

# Delete all variables utilized
del find_junction, get_elevation, get_demand, get_x_coord, get_y_coord

Any advice would be very much appreciated!

Thanks!


texts: Copying data from one table to another automatically using PostgreSQL

texts:
I have a database with two schemas, "public" and "versioned".
In "versioned" I created a copy of the data contained in "public" that I update daily.

Is there an automatism for the data that I add every day on the tables in "versioned" to add them to the tables in "public" automatically on a daily basis?


texts:
I would like to work with selections. For that I start to read the layer from a mxd. I have controlled that the layer count 0 returns the correct name, but the GetCount_management is empty as is the following selection. I have never used the GetCount_management before, it is just a step I added before my selection to check why I do not get data returned. The layer is a shape and definitely not empty. It is the first script that I start in 10.6.1

Stationer = arcpy.mapping.ListLayers(mxd,"",df)[0]
print Stationer.name
print arcpy.GetCount_management(Stationer)


texts:
In QGIS, I'm trying to create a vector grid layer with cells 150m x 150m, yet it only seems to be possible in degrees.

Is there a setting I need to change, or is there an equivalence I can use to convert degrees to meters?


texts: How to add custom Feature attributes to Shapefile using Python?

texts:
I am seeking a way to take an existing Shapefile that has a Feature set of 200 countries. Each country Feature has an attribute of "NAME." My objective is to create a Python script that adds an arbitrary (for now) additional attribute, say, "POPULATION".

Of course I have the OSGeo and GeoDjango modules installed. I'm as far as:

 from osgeo import ogr

    infile = ogr.Open('sample.shp', 1) #'sample.shp' is a pre-existing ESRI shapefile described above
    inlyr = ogr.GetLayerByIndex(0)

Am I missing an OGR function that will allow me to insert Feature attribute fields into an existing Shapefile?


texts: qgis2web: In which files are the zoom buttons?

texts:
Can anyone help me out in which files (in an exported qgis2web project with OpenLayers3) the zoom buttons are specified? Or directly, the geolocation button? I do a project for school with qgis2web, and for the layout and data presentation i had to edit lots of javascript files after the export, so re-exporting the project is not an option, unfortunately.
I did a seperate export, compared the files and managed to add the geolocation function after the export, it does work if i do hit the button. But the button has no mark, and its location in the window isn't lined up with the zoom in and out buttons.

Now i'm unable to find the right file in which i could line up the button better, or exchange it with another symbol/grafic.


texts: Georeferencing old British army map of Palestine?

texts:
I'm using the georeferencer to try and reference a British army map of Palestine. I want it to later be added on Google earth.

I'm pretty sure it should be Palestine grid 1923 the map ends up appearing distorted in all sorts of places around the Middle East but not in Palestine.

What projection information do I need to input in the georeferencer?
I'm currently entering palestine grid when inputing the map into the georeferencer. I'm adding the points according to the info on the legend on the right side of the map. exporting to WGS 84 / Pseudo Mercator


texts:
I have created an ITable but I need to be able to delete it if it already exists. How can I delete an ITable? This is the code that I have so far which checks for an existing ITable.

public static ESRI.ArcGIS.Geodatabase.ITable CreateTable(ESRI.ArcGIS.Geodatabase.IWorkspace2 workspace, System.String tableName, ESRI.ArcGIS.Geodatabase.IFields fields)
    {
        // create the behavior clasid for the featureclass
        ESRI.ArcGIS.esriSystem.UID uid = new ESRI.ArcGIS.esriSystem.UIDClass();

        if (workspace == null) return null; // valid feature workspace not passed in as an argument to the method

        ESRI.ArcGIS.Geodatabase.IFeatureWorkspace featureWorkspace = (ESRI.ArcGIS.Geodatabase.IFeatureWorkspace)workspace; // Explicit Cast
        ESRI.ArcGIS.Geodatabase.ITable table;

        if (workspace.get_NameExists(ESRI.ArcGIS.Geodatabase.esriDatasetType.esriDTTable, tableName))
        {
            //DELETE ITABLE HERE
        }
    }


texts: GDAL warp & GDAL contour inconsistent results

texts:
I am creating 2 output contours based on the same input DEM, one is clipped to a smaller grid feature (green line) while the other is clipped to a larger irregular polygon (blue line).
The problem is there is a .05m gap between the 2 results and the ends of the smaller feature have a bend to them, where the larger feature is straight.
I am using this code to generate the results where the cutline_layer is either the grid feature or large irregular polygon:
    # clip the DEM from cutline layer
    dem = gdal.Warp(raster_out,
                    dem_file,
                    format='GTiff',
                    xRes=2,
                    yRes=2,
                    cutlineDSName=cutline_layer,
                    dstSRS='EPSG:3424',
                    cropToCutline=True)
    dem = None  # remove from memory

    # create the contour
    temp_shapefile_dir = os.path.join('temp', 'contour')
    out_shapefile = os.path.join(temp_shapefile_dir, f"{id_property}.shp")
    os.system(f'gdal_contour -a elevation -i 2 "{raster_out}" "{out_shapefile}"')

Output

Blue line - Contour from large Irregular cutline
Green line - Contour from smaller Grid cutline
Gray dash - Cutline being used for smaller DEM (green line)


texts: Conversion problem from Lambert93 to WGS84 - R

texts:
I am dealing with IGN shapefiles of french "communes" that are in Lambert93 projection.

temp <- readOGR(dsn = chemin)
 temp@proj4string

which gives

CRS arguments:
 +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80
+units=m +no_defs

I would like to focus on the attached dataset and to use the coordinates of some points :

data.temp <- as.data.frame(temp@data)
  coordinates(data.temp) <- ~X_CHF_LIEU+Y_CHF_LIEU

if we look at the coordinates :

> coordinates(data.temp.centr)
      X_CHF_LIEU Y_CHF_LIEU
0           1743      68564
1           7919      63428
2           7747      63253
3           5284      64353
4           9895      68015
5           7931      63483

But I would like to have this in the WGS84 projection system. In order to do so, I first assign the initial CRS to spatial data frame I extracted. And then re-project it.

data.temp@proj4string <- temp@proj4string
  data.temp <- spTransform(data.temp,
                                  CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +"))

But unfortunately, it gives unrelevant coordinates :

> coordinates(data.temp.centr)
      X_CHF_LIEU Y_CHF_LIEU
0      -1.375898  -5.554938
1      -1.335507  -5.584978
2      -1.336522  -5.586134
3      -1.352299  -5.580093
4      -1.324725  -5.555561
5      -1.335450  -5.584630

but obviously, it is not in France anymore.. Any idea why I'm not ending with proper WGS84 coordinates ?


texts: Spatial subsetting SpatialPolygonsDataFrames completely within selecting area using R

texts:
I am trying to spatially subset in R:

the forest stand data (many small polygons) by the
watershed data (one polygon)

My datasets are both in SpatialPolygonDataFrame format. I don't want to intersect or clip the data, because this creates small polygons on the border. Instead, I would like to select polygons completely within or within, such Select Layer By Location tool in ArcGIS. The output, unlike in ArcGIS, should be a new shapefile.

Here is how it is looks like:
forest stands = green,
watershed: blue border,
selected stands completely within: yellow.

However, I cannot find, how to make this in R?

Using sp::over(x = watershed, y = forest) the output is a a data.frame, not a spatial feature. The tools sf::st_within(x, y, sparse = TRUE, prepared = TRUE) seems promissing, but they operates on sf objects, not on spatialPolygonDataFrame.

Is there any quick spatial subset in R solution that I have overlooked?


texts: Selecting max value of a select by location with ModelBuilder and Python?

texts:
I am new to ModelBuilder and Python scripting.

I am trying to create a model which selects each feature within a dataset, uses a select by location to select all features within a 50m radius, and from the selection results, return the maximum value in particular field.

I have setup a model:

And the settings + python script I am using in the Calculate field is:

fc = "Sample_Area_Layer"

field = "nNCC_1"

Cursor = arcpy.SearchCursor(fc)

row = Cursor.next()

MaxValue = 0

CurrentValue = 0

for row in Cursor:

    CurrentValue = row.getValue(field)

    if CurrentValue > MaxValue:

        MaxValue = CurrentValue

    Cursor.next

The results which are being returned are inconsistent. Some have values which are not in the selection, others have values which are not the maximum value of the selected values - so although the Maximum Value may infact be 9 the formula is only picking up a maximum value of 4.


texts: How to join two polygon layers into a resulting layer with overlapping polygons merged together (including attributes)

texts:
I am trying to find a method (in QGIS and/or GRASS) to join two layers of polygons into a resulting layer with merged overlapping polygons. It's hard to explain in words, so I've drawn a simplified graphic

Has anyone come up with a working solution for this need?


texts: Obtaining city model for modification in CityEngine from LiDAR point cloud

texts:
In ArcMap I transformed the point cloud into Multipoints and then into Singlepoints and I obtained this model in ArcScene.
Is it possible to transform these points into 3D objects?
Finally, I would like to get a 3D city model for modification and CityEngine


texts: How to fix this spatial query in GeoDjango?

texts:
I am new to Geodjango. and I am following the steps given in https://docs.djangoproject.com/en/dev/ref/contrib/gis/tutorial/#

and the steps i followed before https://docs.djangoproject.com/en/dev/ref/contrib/gis/tutorial/#spatial-queries were working fine but after that I got error in the

qs = WorldBorder.objects.filter(mpoly__contains=pnt_wkt)

in this Statement.

Error it is showing is

FieldError: Cannot resolve keyword 'mpoly' into field. Choices are: area, fips, geom, id,
iso2, iso3, lat, lon, name, pop2005, region, subregion, un

and if I am executing this line

sm = WorldBorder.objects.get(name='San Marino')

Then Error is

DatabaseError: column world_worldborder.geom does not exist
LINE 1: ...ld_worldborder"."lon", "world_worldborder"."lat", "world_wor...

Kindly Help ..

Thanks


texts: Can ArcGIS Desktop 10 Add-in Use IronPython?

texts:
Coming from the same background as the original poster in this thread, I wonder if anyone has had any success creating an Arc10 Add-In using IronPython.

I have been able to access ArcObjects via IronPython using the info in the linked thread, assorted blog postings, and a few posts on the ESRI forums.  These neophyte attempts have all been run from the command line and I would to be able to gain access through either a toolbar or a button.

One thought that I had was to create the add-in using the .Net 3.5 template, add a python file, and then call a function in the python file from the button's C# file.  It appears though that add-ins are not yet supported in .Net4?

Here is the code from the .cs file which controls button functionality.

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

//I need to then add references to IronPython
using IronPython;
using IronPython.Hosting;
using Microsoft.Scripting;
using Microsoft.Scripting.Hosting;

namespace PyButtonTest
{
    public class ZoomToLayer : ESRI.ArcGIS.Desktop.AddIns.Button
    {
        public ZoomToLayer()
        {
        }

        protected override void OnClick()
        {
            var ipy = Python.CreateRuntime();
            dynamic test = ipy.UseFile("ZoomToLayer.py");
            test.ZoomToActiveLayerInTOC(); //Call this function in the py file

            ArcMap.Application.CurrentTool = null;
        }
        protected override void OnUpdate()
        {
            Enabled = ArcMap.Application != null;
        }
    }

}

As far as I can tell the Add-in format takes care of getting the application hook.  I should then be able to, in the python file, after referencing and importing all of the necessary libraries get the map document with (I think?)

mxDoc = ArcMapUI.IMxDocument.ActiveView

From here I am stumped.  Google has not turned up anything helpful.  Has anyone else tried this?  Got something to work?  Do I need to wait until add-ins are supported in .Net4 to go this route?


texts: Create a perpendicular line of points to another line

texts:
i have two problems over here. QGIS

1- I need to create a perpendicular line to connect multiple points to multiple lines but these lines have to follow a sequence of column.

2- And these perpendicular lines can't be created of differents "id's".

In the image above, you can see the result (In red) that i want. The (Red Line) has to start on the point and connect to the "line.shp" (following column "MDE" from lowest to highest) creating a perpendicular line.


texts: Expanding raster within mask by interpolating it with boundaries set to 0

texts:
I'm having a hard time trying to expand a raster within a mask. Here are some details about my case:

I'm working with glacier ice thickness rasters, which happen to be smaller than the real size of the glaciers, specified by a vector polygon which I'm using as a mask. My goal is to expand the ice thickness raster within the polygon by interpolating the raster with the borders of the polygon which should be 0 (meaning no ice). I've tried to address this problem in multiple ways, and I have ended up converting the raster to points and the borders to points with value 0 and merging both layers. Here is what I have so far tested:

When I try to interpolate the points using different interpolation tools I'm getting this error message: 'NoneType' object has no attribute 'crs' See log for more details

I have also tried doing the same but directly with a raster file (with the ice thickness raster added surrounded by 0 values), and trying to interpolate it with the glacier polygon mask.

I've been using GRASS' r.surf.idw tool, but I can't get it to work. Despite having the raster file and the polygon file with the same CRS, I see a warning label on top of my screen saying Warning: The projection of the chosen layer is not the same as canvas projection! The selected extent might not be what was intended.

The results always look like this:

Anyone has any ideas on how to tackle this issue?


texts:
I'm a database administrator who's working with a new client who'd like to use a "search" or "filter" application with GIS. For example, he would like to search within his area of properties that are over 30 acres, have this elevation, etc.

Is there a way to search in that system or to create a search/filter system?


texts: Algorithm for Cost Connectivity tool of ArcGIS Spatial Analyst?

texts:
For the Cost Connectivity tool in the Spatial Analyst toolbox, what is the algorithm used to compute the minimum spanning tree?
The ArcGIS website includes a general description of how the cost connectivity algorithm works, but it does not specify what minimum spanning tree algorithm is used.
I am asking because I am a PhD student who may use the cost connectivity tool for my research, and it is important that I understand how the tool works.
Edit to clarify based on responses so far: Please note I am specifically interested in the minimum spanning tree algorithm (e.g., Dijkstra). From the Esri website, I understand generally how the Cost Connectivity tool otherwise works.


texts: Georeference multiple raster layers of same area at once or identically

texts:
I have many raster layers that I need to georeference to my existing vector layers.  Data types vary between 8 bit unsigned, 32 bit float, and bit 64 float. All have the exact same extent, cell size, and from and to areas. Is there a way to georeference all the layers at once so they all have the same RMS error? The area is complicated as we do not have roads etc. for exact points of reference.

I am using QGIS.


texts: Avoiding setting env.workspace and other paths/settings too often in ArcPy?

texts:
I have some scripts that are constantly setting env.workspace of arcpy. Usually, I have a caller script(script1.py) that may call some other script (script2.py). I've created a file globals.py with the paths that never change:

import arcpy.env as env

rootPath="S:/workspace/MSc_git/ArcGis/"
layersPath=rootPath+'Layers/'
# Set the workspace environment to local file geodatabase
env.workspace=rootPath+'estagio_db.gdb'

In script1.py I have something like this:

import arcpy
import globals as g
from arcpy import da
arcpy.AddMessage("Importing arcpy outside ArcMap took %s seconds." % (time.clock() - t_start))
from cProfile import run

def main():
  ....

And I need to use env.workspace and sometimes the layerspath inscript2.py
How can I do this (should I import gloabl.py as well?)? should I pass it by argument when calling a function from script2.py in script1.py or is there any other way to accomplish this?

I know we should avoid global vars, but in this case, and since I need to centralize its definition I don't know how to do it well...


texts: Are there scenarios where the JSON payload from an AGOL webhook would have multiple collections in the array?

texts:
I have a webhook in AGOL that posts a JSON payload to an external system (Integromat).

The webhook is triggered when features are edited in a feature layer (Change Types: FeaturesEdited).

This is what the JSON payload looks like:
[
   {
      "name":"Webhook 1",
      "layerId":0,
      "orgId":"Ak",
      "serviceName":"Work",
      "lastUpdatedTime":1597279482838,
      "changesUrl":"https://services.arcgis.com/Ak/arcgis/rest/services/Work/FeatureServer/extractChanges?serverGens=[4,7]&async=true&returnUpdates=false&returnDeletes=false&returnAttachments=false",
      "events":[
         "FeaturesCreated"
      ]
   }
]

Question:
I've noticed that the JSON collection in the payload is wrapped in an array.
Are there possible scenarios where multiple collections would be included in the array?

Or can I safely assume that there will only ever be a single collection in the array and always point to the first colleciton?


texts: gdal_merge runs indefinitely, crashes system, BSOD

texts:
I have two small rasters. Same resolution, SRID, type, format. When I use gdal_merge, the console keeps running the command ad infinitum.
import time, os

ot = "UInt16"
of = "GTiff"
nodata = "255"
ps = "5 5"
output = r"D:\Tasks_D\merge_minpnt_bris_yate.tif"
input1 = r"D:\Tasks_D\Raster_dev\gdal_sql_minpnt_rec_bris_5m.tif"
input2 = r"D:\Tasks_D\Raster_dev\gdal_sql_minpnt_rec_yate_5m.tif"

merge = fr'gdal_merge.py -n {nodata} -a_nodata {nodata} -ot {ot} -o "{output}" -of {of} -ps {ps} "{input1}" "{input2}"'
print(merge)
os.system(merge)

Here is what my console looks like after I run, note that I have to kill the script otherwise it keeps repeating the console command until I get memory issues and even Blue Screen of Death (BSOD):

Same issue occurs even if I go the gdal_buildvt route. Either way, I've tried to let this run and nothing happens, just repeating code. This should only take a second or two. I've tried to run this on two other tifs, sometime overlapping tifs, sometimes ones adjacent to each other. No luck. Any ideas on what's happening?
FYI I have a powerful PC with 64GB RAM.


texts: Is it possible to get the EPSG value from an OSR SpatialReference class using the OGR Python API?

texts:
When reading a layer from a OGR PostGIS connection I can get the SpatialReference of the layer, but is it possible to get the EPSG value? Is there any documentation on this?

For example:

lyr = conn.GetLayerByName(tbl) # Where conn is OGR PG connection
srs = ly.GetSpatialRef()
print srs

Returns:

PROJCS["OSGB 1936 / British National Grid",
GEOGCS["OSGB 1936",
    DATUM["OSGB_1936",
        SPHEROID["Airy 1830",6377563.396,299.3249646,
            AUTHORITY["EPSG","7001"]],
        AUTHORITY["EPSG","6277"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4277"]],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",49],
PARAMETER["central_meridian",-2],
PARAMETER["scale_factor",0.9996012717],
PARAMETER["false_easting",400000],
PARAMETER["false_northing",-100000],
AUTHORITY["EPSG","27700"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]

So how do I get the EPSG value for the projection? E.g.:

srs.GetEPSG()
print srs
27700

I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG'.


texts: Bypassing 10 character limit of field name in shapefiles?

texts:
I'm exporting geometry with attached text attributes from oracle database to esri shapefile format (.shp) with Java and Geotools library.

Attribute columns in our database have names with more than 10 characters, and Geotools forces to truncate them. I understand that it's due to .shp or .dbf files specification.

I can workaround this, by creating some simple txt file with "shrtname"="The full and long name", but obviously it won't be understand and imported by any other software than ours.

Is there official way to do the mapping from short field names to long full-text names?

For example, xml file next to all other .shp .dbf .shx files.


texts: Getting values of last row in table with ArcPy/SearchCursor?

texts:
I would like to base my query on the last row of a table?  I know the basic syntax/functionality for my search cursor, how do I query only the last row;

Below is an example of my search criteria,  this data comes from a log of success/fail for a script that runs, I will base other decisions off of the success or fail later in the application.

How do I add to my existing code the function to look in the last row, to getValue of timeand success or fail

I am using ArcGIS 10.2.

cursor = arcpy.SearchCursor(aTable,"""Time > DATEADD(minute, -2,  GETDATE())""")
for row in cursor:
    print(row.getValue("Time"))
if row.getValue("Time") < datetime.datetime.now().strftime('%m/%d/%Y %H:%M:%S'):
    print "true"
else:
    print "false"

This is later implemented in the script and creates the success or fail values as yes and norespectively.

rows = arcpy.InsertCursor(aTable)
for x in xrange(0, outputobjects):
    row = rows.newRow()
    row.setValue("Success", 'YES')
    row.setValue("Fail", 'NO')
    row.setValue("Time", Start)
    rows.insertRow(row)


texts: Rendering problems for some polygons in a Leaflet map with a GeoJSON tile provider

texts:
New to mapping.

Trying to create a heat map with data that updates every hour for US Counties with leaflet.js similar to http://leaflet.cloudmade.com/examples/choropleth-example.html.  But, adding all the US Counties in 1 shot makes map jerky in IE.

So found this example of using GeoJSON tiles with leaflet.js, http://bl.ocks.org/3725681

Only problem is that many tiles are not rendered completely.  I think some of the z/x/y.json files are missing coordinates, but not sure.  Can someone point me to the right direction of the problem here.


texts: Edit KML with the mouse in QGIS: is it possible?

texts:
In google maps, and more recently in Google Maps Engine Lite, you can create a Track/LineString by clicking with the mouse along a route, delete some intermediate point, drag the midpoint of a segment to decompose it in two subsegments, and so on.

This is pretty much what the "Nodes" tool can make in drawing programs like CorelDraw, Inkscape and Illustrator.

Is there a way to have such functionality in QGIS?

My typical use-case would be to open a KML file with pre-existing tracks, modify those tracks by adding, moving and deleting points (and ideally by joining and splitting tracks), and save it back to the original file.

I am very new to QGIS, and if I can make this it would be the perfect software to my needs.


texts: How to handle PostgreSql 1:m relationship data by using “Data-Driven Input Mask” QGIS plugin?

texts:
I have 1:m relationship tables in my PostgreSql. I like to show the main table attributes as well as the relational table attributes (like grid rows) in QGIS attribute-edit-form.

What are the entries i have to add to dd_* tables to get this?