texts: Unable to change name of feature class using ArcGIS Desktop?

texts: I am in ArcView 9.3.  I have created a personal gdb and a Feature Class.  I was able to change the name of the gdb, but why CAN'T I change the name for the Feature Class (FC)?


texts:
In PyQGIS Cookbook, says, Any background task (regardless of how it is created) must NEVER use any QObject that lives on the main thread, such as accessing QgsVectorLayer, QgsProject or perform any GUI based operations like creating new widgets or interacting with existing widgets. But I have to use QObject to give feedback to User. For example, I have a Qlabel that keeps the user up to date on how many layers have been processed. If I change the text of label in thread, Qgis crashes. How can I change the label's text during processing layers. I use this approach:
btn_run.clicked.connect(task)
def task():
    t1 = threading.Thread(target=compilePage, args=("compilePageThread",))
    t1.start()
def compilePage(threadName):
    label.setText('Please Wait.')
    for layer in layers:
        #do some work on layer and then change label text to give feedback
        label.setText('{} has been processed.'.format(layer.name()))


texts: How to import Vietnamese characters with MapInfo Pro?

texts:
How can I import shapefiles, Charset "UTF-8" with Vietnamese characters (f.e. "Đường nhựa") correctly?

With UTF-8 I get something like "___ng nh_a"!

Should I use another charset or is it just not possible with MapInfo? QGIS works fine with the same shapefile.


texts: Determining the slope of a DEM (SRTM-data) in QGIS

texts:
I'm new to QGIS and I have to calculate the slope of a DEM (I use srtm-data). I already detected the geoalgorithms toolbox and I also produced a raster layer with the slope using the tool "slope, aspect, curvature" (One can find it in Domain specific algorithms - Terrain analysis and geomorphometry). I used the 2. Degree Polynom (Zevenbergen & Thorne), which is the default. A visual check shows that the general location of slopes seems to fit (I compared it to Google Earth).

Now my problem is: I have no clue about the unit of the output. The pixels at the steepest slopes have values around 0.3.... I have never been to that area, but it is located in Southern Germany, in an area that is usually quite hilly. So I strongly doubt that the slope is only 0.3° in the steepest area. I have the assumption that it is just 30°, but I would like to check somewhere to be absolutely sure how to interpret the result (since it could also be percentage values). I find a lot of information about what the tool does, but never any site where the units of the output are clearly stated.

Does anybody know how I can find out the unit of this output?


texts: Importing KML Google Maps files with custom data tables to QGIS?

texts:
I am trying to visualize the complete attribute table of an exported Google Maps KML file with a customized data table, but QGIS 2.4.0 won't parse the added information correctly.

This is what I can read in the data in text format (1 line of the table):

 <name>Piante aliene Filipponi-Raniolo, Naviglio Pavese (zona 10)  06/06/2015</name>
    <description><![CDATA[]]></description>
    <Folder>
        <name>Livello senza titolo</name>
        <Placemark>
            <name>Marker 1</name>
            <description><![CDATA[Acer negundo: 1
Ailantus altissima: 1
Robinia pseudoacacia: 1
Parthenocissus quinquefolia: 1
Euforbia: 1
Acetosella: 1
Superficie verde/tot: 1]]></description>

            <ExtendedData>
                <Data name='Acer negundo'>
                    <value>1.0</value>
                </Data>
                <Data name='Ailantus altissima'>
                    <value>1.0</value>
                </Data>
                <Data name='Robinia pseudoacacia'>
                    <value>1.0</value>
                </Data>
                <Data name='Parthenocissus quinquefolia'>
                    <value>1.0</value>
                </Data>
                <Data name='Euforbia'>
                    <value>1.0</value>
                </Data>
                <Data name='Acetosella'>
                    <value>1.0</value>
                </Data>
                <Data name='Superficie verde/tot'>
                    <value>1.0</value>
                </Data>
            </ExtendedData>
            <Point>
                <coordinates>9.1735411,45.4334849,0.0</coordinates>
            </Point>
        </Placemark>

Right now, if I try to open the KML in QGIS directly, I only get the two default columns (name, description) with description containing all information.


texts: Feature class is not imported in ArcSDE in ArcCatalog

texts:
I installed ArcEditor and ArcSDE. In my computer ArcEditor 9.3.1, ArcGIS Server 9.3.1, ArcSDE 9.3.1 and Oracle 11g are installed. I did all things by the Guide of ESRI: ArcGIS Server, Configuration and setting for Oracle. All the procedures passed normally.
But when I import feature class in ArcSDE in ArcCatalog, the process starts but does not end. I tried to copy the class. Window of primary parameters opens. After clicking OK a window of copying opens. But it is blank and hangs.


texts: Temporary layer has no attribute '__getitem__' with pyqgis

texts:
I used the polygonize tool in order to polygonize a raster layer and it worked well, but now, I would like to save this new vector layer as a temporary layer.

My code from various posts is the following :

grid = QgsVectorLayer("Polygon?crs=epsg:4326",
                      "Grid_"+str(self.id),
                      "memory")

#grid = self.pabs+"/1_Grids/Grid_"+str(self.id)+".shp" #non temporary layer initially used

layer = self.pabs+"/0_Inputs/LandUse/lu_{}.tif".format(self.id)
processing.runalg('gdalogr:polygonize',layer,"lu_majorit",grid)

When I execute the code in the python console of Qgis I get:

'QgsVectorLayer' object has no attribute '__ getitem__'

I know that you sometimes have this problem with python class but here I don't see where the problem comes from.

I was wondering if it is linked to the fact I create an empty layer ?

I am using QGIS 2.18.10


texts: create polygons and fishnets of specific size and orientation and buffers

texts:
I want to study 6 areas in the field and in order to do so I want to divide each one of them into smaller areas/parts.

What I would like to do is find an average square or rectangle area surface (more particularly I need both figures) and then add 4 buffers of 50m each (meaning a total of 200m buffer zone), in each of my study areas. All surfaces could be made using fishnet of 50mx50m too. I guess I can add data in each polygon / fishnet cell, by editing the attribute table in each polygon after it is created.
I am not a GIS expert, so indicate the steps in order for me to reproduce each procedure (creating shapefiles, editing vector, making the buffer etc).

I have georeferenced google earth images of my areas using QGIS 2.4. I have done a little search in this site for possible answers, but quite many of these solutions do not work in the QGIS I use; either errors occurred or do not suit my case.

I am adding some pictures to give a visualisation example of what I would like to do, made by using a presentation (powerpoint) program. I am giving 2 examples for each case, rectangle or square. With blue is my main area and with orange the buffer zone(s). In case of fishnet, every eye should be 50mx50m. In the other case (if fishnet is not possible) there could be polygons with the same snapping point and a 50m distance in between each polygon's side with its closest one. Both  surfaces, rectangle or square, enlosed by blue line should cover the same area in sqm.

below you can see the examples:

rectangle with fishnet
rectangle with polygons (if fishnet option is not applicable)
square with fishnet
square with polygons (if fishnet option is not applicable)


texts: arcpy.AddMessage does not produce newlines in ArcGIS Pro

texts:
I would like to produce additional newlines when using arcpy.AddMessage .  I am using arcpy code assigned to a geoprocessing tool and want to have a newline between some messages. This doesn't happen.

import arcpy

line1 = "This is line 1 \n"
line3 = "This is line 3" + "\n"
line6 = "\nThis is line 6"

arcpy.AddMessage(line1)
arcpy.AddMessage(line3)
arcpy.AddMessage(line6)

Produces

This is line 1
This is line 3
This is line 6
Completed script Script...

I've seen the answer to Using Python new line (\n) with ArcGIS geoprocessing progress dialog?, and the answers there indicate that the above should work.  I note that the OP is using ArcGIS for Desktop. I am using ArcGIS Pro 2.4.1.


texts: Longitudinal and Latitudinal coordinates of interstate exits

texts: I need to find the coordinates of all the interstate exits for some research.  I'm not looking to code up a program that will do it for me, but I was hoping that someone had already done it and would be willing to share that information with me.  Any help would be much appreciated.


texts:
I have a csv file with the positions (x and y) outlining multiple complex polygons. Each point is numbered in a separate field and each polygon is numbered in a separate field. I can load those points into QGIS easily but I cant convert them into distinct polygons. I know there are multiple similar questions but nothing works - convex hull crashes and minimum bounding geometry is too basic. I havent found a distinct package either (im using QGIS 3.14.16).
I dont know. I can I transform a list of points into multiple polygons ? Should be easy but its not.


texts: Which ESRI defined geographic transformation do I need to use to convert my " NAD83 Northern Virginia" data into WGS 1984 Web Mercator?

texts:
My title pretty much tells it all..

I have a number of datasets in which I want to batch project from NAD83 Northern Virginia Datum into the correct Web Mercator projection..

When I do a straight re-projection I get an error message.. this tells me that I need to use a transformation of some kind.. I think..

So, here in lies my problem.. as there is not a web mercator projection transformation listed in the geographic transformation I was thinking I would just use a generic one like...

NAD_1983_To_WGS_1984_2 (there are 8 different ones)
WGS_1984_Major_Auxillary_Sphere_To_WGS_1984

But I am not sure that is correct.. Which is the correct "geographic transformation"?


texts: Adding file name to attribute field using Calculate Field?

texts:
Does anybody know how to add the shapefile name to a field in the file's attribute table?

I've found a description:

  'You have to use inline variable in CalculateField like %Name% in
  expresssion so that whatever the name of the file is, the string will
  be inserted based on Name of the dbase file !'

I've tried to CalculateField tool with %Name%, but I got message:

Executing: CalculateField "2011-11-6 15_41_point" DBFName %Name% PYTHON #
Start Time: Thu Dec 08 23:34:23 2011
WARNING 000405: No records within table
Succeeded at Thu Dec 08 23:34:24 2011 (Elapsed Time: 1.00 seconds)

Where is my mistake?


texts: Creating parallel lines in MULTILINESTRING table using Postgis?

texts:
Possible Duplicate:
  How to create one-sided buffers or parallel lines in PostGIS?

I need to make parallel lines on a MULTILINESTRING table, similar to st_buffer, but that the result come be lines, not polygons. I saw some posts about the st_offsetcurve function, and I thought this is the solution to my problem, but I didn't see this function on Postgis 1.5! Sorry for my question, I'm a new Postgis user!


texts: Rasterizing arbitrary geometry to GeoTIFF in Python?


texts: Models on QGIS 3 are stored following a XML schema. Do you know if it is possible to embed metadata, may be using Dulin Core or ISO? I have noted that the modeler's help take in account subjects like: author, version, documentation and so on. But I am not sure how to include metadata information on .model3 files.


texts: Error in merge function "Erreur de géotraitement GEOS: la différence a échoué"

texts:
I am trying to merge three large shapefiles. After cleaning all the inconsistencies and artefacts, the merging of file 1 and 2 ends well, yet with the following message (sorry to be in French):

  L'entité n'a pas pu être écrite
  suroutput_126be82a_bc81_4c47_908a_c1c95da91787

  Résultats: {'OUTPUT': 'output_126be82a_bc81_4c47_908a_c1c95da91787'}

I don't think it is a big problem.

But then, when I try to merge the 3rd file, the same error pops up and the algorithm stops with the following error:

  Erreur de géotraitement GEOS: la différence a échoué.

and no more details are given.


texts: Create feature class with ArcCatalog in Oracle Spatial without ArcSDE error

texts: I need to create a table/feature class with ArcCatalog 10.2 directly to a Oracle Spatial 11.2 database without ArcSDE, Is it possible? I can execute the create new feature class in ArcCatalog, but after I finish all steps of the wizard (I created a connection to the database with ArcCatalog (Add Database Connection), and clicked right button over that connection -> New -> Feature class) I have the following error: "Failed to create feature class. Undelying DBMS error [ORA-24344]" Could be a database problem? Perhaps something about privileges?


texts: Cannot add Plugins from Python: "Plugin is Broken"

texts:
I'm using the latest version of QGIS on a Windows 7 (32 bits), and when I try to add a plugin from the option "Fetch Python Plugins", it simply cancels the installation and I can see the following error:

Couldn't load plugin 'xytools' from ['C:/PROGRA~1/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/ProprietÃ¡rio/.qgis//python', 'C:/Users/ProprietÃ¡rio/.qgis//python/plugins', 'C:/PROGRA~1/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\Program Files\ArcGIS\bin', 'C:\PROGRA~1\Quantum GIS Lisboa\bin\python27.zip', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\DLLs', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\plat-win', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\qgis\bin', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\Quantum GIS Lisboa\apps\qgis\python\plugins\fTools\tools']

The plugin that was not successfully installed appears in red colour saying "Invalid" - and in the description field: The plugin is Broken.
It happens with all plugins from the QGIS Plugin Installer.
Does anyone knows how can I fix this problem?

Thank you in advance for any help you can provide.


texts:
Possible Duplicate:
  How can I remove a feature layer in an Arcpy Script

This question could seem too stupid, but i'm struggling with it, and it cannot work.

I have an mxd file containing raster layers, and I would like to remove a layer.  I don't know why it does not remove, it doesn't pop up any errors.  I'm trying this code:

mxd = arcpy.mapping.MapDocument("photo.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]

for i in arcpy.mapping.ListLayers(mxd , "TerraColor_SanFrancisco_US_15m.tif"):
    arcpy.mapping.RemoveLayer(df , i)

arcpy.RefreshActiveView()

Any help please.


texts:
I have a custom tool to add feature points. The attribute values are provided filling a custom dialog.
The code was writted following this topic.

At the end of the procedure, outside the tool, if I select a different layer and I start editing with toggle editing of QGIS (without selecting add point/polygon feature tool), I continue editing with the custom tool as if it didn't end the procedure.

How I can stop the editing session on the custom tool?
This is the code:

class casePicker(QgsMapTool):
    afterClick = pyqtSignal()

    def __init__(self, dlg, psrid, iface, tt, lab, lyr):
        QgsMapTool.__init__(self, iface.mapCanvas())
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
        self.dlg = dlg
        self.psrid = psrid
        self.tt = tt
        self.lab = lab
        self.lyr = lyr
        self.funcs = qvfuncs.VetEpiGISFuncs()
        self.canvas.setCursor(mutato)

    def canvasPressEvent(self, event):
        pass

    def canvasMoveEvent(self, event):
        pass

    def canvasReleaseEvent(self, event):
        x = event.pos().x()
        y = event.pos().y()
        pt = self.canvas.getCoordinateTransform().toMapCoordinates(x, y)
        x = pt.x()
        y = pt.y()

        if self.dlg.exec_() == QDialog.Accepted:
            QApplication.setOverrideCursor(Qt.WaitCursor)
            x = float(self.dlg.lineEdit.text())
            y = float(self.dlg.lineEdit_2.text())
            self.addFeat(x, y)
            self.tt.setChecked(False)
            QApplication.restoreOverrideCursor()

    def addFeat(self, x, y):
        self.lyr.startEditing()
        feat = QgsFeature()
        feat = self.funcs.outattrPrep(self.dlg, self.lyr)
        pnt = QgsGeometry.fromPointXY(QgsPointXY(x,y))
        if self.lyr.geometryType() == QgsWkbTypes.PointGeometry:
            feat.setGeometry(pnt)
        else:
            pass
        feat.setValid(True)
        self.lyr.addFeature(feat)
        self.lyr.commitChanges()
        self.lyr.updateExtents()

    def activate(self):
        pass

    def deactivate(self):
        pass

    def isZoomTool(self):
        return False

    def isTransient(self):
        return False

    def isEditTool(self):
        return True


texts:
Has anyone had a similar "extreme slow" experience using "Select by location" (Toolbox) and "Spatial query" tool?

I've used both for selecting some points inside a polygon. The data was in SHPs and the point-dataset has 80.000 points (with a lot of attributes) - the poly-dataset is small.

The result was "within" ~14.000 points. QGIS 2.8 and 2.13 need around 10min for the task. I did the same now with ArcGIS 10.3 and it took 4-5 seconds.

Does someone know a better (faster) solution within QGIS ?


texts: Stop queued click event if another event is fired

texts:
I'm looking for a way to suppress a map click event if some other event has already happened. The application I have created has two primary functions:

Click on marker (park points) --> Brings up information about the park you clicked on
OR
Click on the map (outside of a marker) --> Clears the map of all information and graphics and runs a query for parks within the area and spits out a list

Here's what I've got:

var hover_toggle = false;

...

// Map event handlers
click_event = map.on("click", function(evt){
  if (hover_toggle == false) {
     var point = evt.mapPoint;
     var graphic = evt.graphic;
     parkQuery(point, graphic);
  }
});

// Marker event handlers
parks.on("mouse-over", function(evt){
    map.setMapCursor("pointer");
    hover_toggle = true;
});
parks.on("mouse-out", function(evt){
    map.setMapCursor("default");
    hover_toggle = false; });
parks.on("click", function(evt){
    var park_point = evt.graphic;
    parkZoom(park_point)
});

...

This has it so that whenever a user is hovered over a marker, the hover_toggle makes it so the map.on("click") event is suppressed. However, this doesn't always work... sometimes the marker event fires (it seems to always happen first), and then the map event fires and clears the information the user brings up.

Any ideas on how to suppress or disconnect/connect the map.on("click") event if the marker event is fired at the same time?


texts:
I use WAMP in my Windows, and I want to use phpmapscript for programming. I've copied the phpmapscript.dll file to C: \ wamp64 \ bin \ php \ php5.6.25 \ ext and I added the "extension = php_mapscript.dll" in the php.ini file.
In the wamp extension menu, this plugin is checked and added to this list

but when I use the -phpinfo- command , I do not see this plugin.

I have a lots projects inside WAMP and I can not use ms4w on my system.

Can you guide me?


texts: Trace Downstream task: which ArcObject interface and method?

texts: Does anyone know what ArcObjects interface and method the Trace Downstream trace task on Arcmap's Utility Network Analyst toolbar calls?


texts: How to create a crisp topographical/terrain map with ggplot2?

texts:
I am trying to produce a high quality vector based map (in greyshades) in order to report locations.

I tried this by getting maps from GADM and SRTM to make a raster as a base for plotting in ggplot2. However, the resulting data frame is too large for plotting. Quesion 1: How to simplify the data frame for ggplot2 without while still obtaining a high quality resolution map?

Here is what I tried for south New Zealand:

library(dplyr)
library(ggplot2)
library(raster)
library(rasterVis)
library(scales)
library(rgeos)

nz1 <- getData('GADM', country='NZL', level=1)
nz1 <- subset(nz1,NAME_1 %in% c("Southland","Otago","West Coast"))

nz1c <- gCentroid(nz1) %>% coordinates()

dem1 <- getData("SRTM",lat=nz1c[2],lon=nz1c[1],path=datadir)
dem2 <- getData("SRTM", lat = -45.516667, lon = 168.566667,path=datadir) # City of Athol coordinates
dem3 <- getData("SRTM",lat = -45.866667, lon = 170.5,path=datadir) # City of Dunedin coordinates
dem4 <- getData("SRTM",lat = -44.383333, lon = 168.716667,path=datadir) # Mount Aspiring GPS data -44.383333, 168.716667

dem <- merge(dem1,dem2,dem3,dem4)

dem <- crop(dem,nz1,filename=file.path(datadir,"dem_nz1.tif")

dem.p  <-  rasterToPoints(dem)
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

p1 <- ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt))

The data frame has about 25 Million observations Which produces something like this (~110MB PDF):

This looks already very promising but it is way too heavy.

How could I possibly reduce the data frame or use any other approach to obtain a nice high scalable vector image for my ggplot2 maps?

I just started working with maps in R and basically I just want to produce nice terrain/topographic maps, and I don't want to use ggmaps or similar for this.


texts: Reprojecting: Different results in OGR and PROJ.4 for EPSG24373 (India)

texts:
I need to convert EPSG4326 lon/lat tuples to EPSG24373 (Kalianpur 1880 / India zone III).  However, when using the Python bindings of OGR and PROJ.4 (pyproj), I get different results.

I wrote the following code:

LON, LAT = 76.661217, 15.183595

import numpy as np
from osgeo import ogr, osr

src_epsg = 4326
dst_epsg = 24373

source = osr.SpatialReference()
source.ImportFromEPSG(src_epsg)

target = osr.SpatialReference()
target.ImportFromEPSG(dst_epsg)

trans = osr.CoordinateTransformation(source, target)

point = ogr.CreateGeometryFromWkt("POINT ({} {})".format(LON,
                                                         LAT))
point.Transform(trans)
print("OGR result: {}".format(point.ExportToWkt()))

from pyproj import Proj, transform
proj = Proj(init="epsg:{}".format(dst_epsg))
x1, y1 = proj(LON, LAT, radians=False)

print("pyproj result: {}".format((x1, y1)))

which yields the results

OGR result: POINT (2607326.263296442106366 542082.358424285310321)
pyproj result: (2384135.3043296584, 495679.3120856479)

From my expectation, the results should be identical, but they aren't.  What am I missing here?


texts: Unable to get soil data for coordinates using Soilgrids REST API

texts:
I am trying to extract soil data from https://soilgrids.org/ using their very new Soilgrids REST API. They have a special website (http://rest.isric.org/soilgrids/v2.0/docs) where I am able to enter the coordinates and get the soil data that I require.

But for some reason, I am not able to get data for many coordinates that I tried.
For example, I can get data for (lat=39.1813855,lon=-77.1827295) but not for (lat=39.1825835,lon=-77.1828465).

Why can't I get data for the last point (lat=39.1825835,lon=-77.1828465)??


texts: Google Earth Engine Reducing Landsat Global Inland Waters Raster to Vector

texts:
I'm trying to create a vector outline of the water areas of the Landsat Global Inland Water raster in Earth Engine. I get the error "Line 3: water.addBands is not a function" from the following code:
var water = ee.ImageCollection("GLCF/GLS_WATER");

var vectors = water.addBands(water).reduceToVectors({
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'regions'
});


texts: Symbolise polygon outline where it doesn't coincide with other feature using geometry generator in QGIS?

texts:
I would like to symbolise a polygon without converting them into a line. I have a polygon feature (hatched and bounded by red line) which coincides on top and bottom of the coastline (blue line). I would like to symbolise the outline of a polygon only on two sides and want to drop out the red colour on top and bottom.

I would like to see the end result like this.


texts: How to save images from Google Earth Engine using Python?

texts:
I would like to export images from image collection using Python. This is what I am doing:
cover = ee.ImageCollection('MODIS/006/MOD11A2').filter(ee.Filter.date('2018-01-01', '2018-05-01'))\
  .select('LST_Day_1km')

Then I clipped the image over a geometry of interest from a shapefile
import geopandas as gdp
myShp = gdp.read_file(myGdp)
def returnCityBoundary(gdf):
  bounds = gdf.total_bounds
  lon_point_list = [bounds[0], bounds[0], bounds[2], bounds[2]]
  lat_point_list = [bounds[1], bounds[3], bounds[3], bounds[1]]
  polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
  crs = {'init': 'epsg:4326'}
  polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])
  geometry = ee.Geometry.Polygon(
          [[[bounds[0], bounds[1]],
            [bounds[0], bounds[3]],
            [bounds[2], bounds[3]],
            [bounds[2], bounds[1]]]])
  return geometry

  geometry = returnCityBoundary(gdf)

In this way I can clip the images over the region of interest.
cover_clipped = cover.map(lambda image : image.clip(geometry)) ## get images only for the region (clip to geometry)

I can also show the images as described here. This an example from Paris:

How should I do to save ImageCollection as images? I am trying to do the following but it only opens the drive folder without saving anything inside.
task = ee.batch.Export.image.toDrive(**{
    'image': cover_clipped,
    'description': 'imageToDriveExample',
    'folder':'Example_folder',
    'scale': 100
})
task.start()


texts: Where to Start: .Net(Xamarin) map viewer with cached map and overlay ability

texts:
I am new to the whole GIS topic. Two days ago, I was asked to develop a map viewer.

I know that

it should be a Xamarin Application
it should display a local cached map (a vector map)
I should be able to draw shapes on it at specific coordinates

The last two days I browsed the Internet to find some info on where to start. I found this Stackexchange, and many other free and licensed developer kits. I think openlayers would suffice my needs, but it's mainly for web applications which is a viable option, if I cannot find something good for .Net.

I found ArcObjects but it is not free and I am not sure if this is what I am looking for. It seems that it is used mainly for ArcGIS Plugins.

Also I found Mapbox SDK, but it only supports iOS and Android.

Has anyone experience with dotspatial?

I hope you guys can point me to the right SDK for my needs.

Edit:

An SDK was not the right thing. I am looking for a good GIS Framework. For displaying maps and drawing shapes onto it.


texts: Making the results from function a string

texts:
I am trying to get the features count from feature, then put it in a dictionary that I can call later:

ct_dict = dict()
...
fipsList = ['13245','06007']
for fip in fipsList:
    mergedJoinCentroids = C:/Workspace/Sandbox/MapChangeProject/selects.gdb/mergeJoinCentroids_{}".format(fip)
    property_ct = str(arcpy.GetCount_management(mergedJoinCentroids))
    ct_dict.setdefault(fip, property_ct)

I tried to make a simple function that I could then easily import into another script:

def ct_dict():
   return ct_dict

This should return a dictionary that looks like this {'13245': '72456', '06007': '93557'}.However, I run the other script with the ct_dict() function, I get an error: AttributeError: 'function' object has no attribute 'items':

import <script name>
from <script name> import ct_dict
...
counts = ct_dict()
...
for k, v in ct_dict.items():
    if fip == v:
       sheet1.write(r, propertyCol, v)

So, basically, it is just returning something like this <function FEMA_WMS_3.ct_dict()>, not a dictionary. How do I resolve this?


texts: Setting polygon style with setStyle() function in Leaflet does not work from select menu

texts:
I think this question will be easy for you
But I don't know where is the problem

I am trying to change the style of the rectangle using the select menu but it does not happen

In contrast, when using the buttons, the change is done

code in "codepen"
https://codepen.io/A_Abdelwahed/pen/yLYoPqj

    <button id="style2" style="position: absolute;z-index: 999;" >style2</button>

    <button id="style3" style="position: absolute;z-index: 999;" >style3</button>

    <button id="style4" style="position: absolute;z-index: 999;" >style4</button>

    <select id="selectDataType" style="position: absolute;z-index: 999;" ></select>


texts:
I would like to translate a vbsript expression into a python expression. But it doesn't work!
Here is the code in vbscript:

Function FindLabel ( [CONTENU], [NOM_OBJET_USUEL], [ALTITUDE_FIL_EAU_M] )
If ( [CONTENU] = "Eaux usées") then
    FindLabel = "<CLR red='168' >" & [NOM_OBJET_USUEL] & vbnewline & [ALTITUDE_FIL_EAU_M] & "</CLR>"
ElseIf ([CONTENU] = "Eaux pluviales") then
    FindLabel = "<CLR green='77' blue='168' >" &[NOM_OBJET_USUEL] & vbnewline & [ALTITUDE_FIL_EAU_M] & "</CLR>"
ElseIf ([CONTENU] = "Eaux mélangées") then
    FindLabel = "<CLR red='79' green='41' >" & [NOM_OBJET_USUEL] & vbnewline & [ALTITUDE_FIL_EAU_M] & "</CLR>"
ElseIf ([CONTENU] = "Eaux de drainage") then
    FindLabel = "<CLR green='168' blue='132' >" & [ALTITUDE_FIL_EAU_M] & "</CLR>"
End if
End Function

My translation:

def FindLabel  ([CONTENU], [NOM_OBJET_USUEL], [ALTITUDE_FIL_EAU_M] ):
 if [CONTENU] == "Eaux mélangées":
    return  "<CLR red='79'  >" + [NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] + "</CLR>"
  elif [CONTENU] == "Eaux usées":
    return  "<CLR red='168' >" +[NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] +"</CLR>"
  elif  [CONTENU] == "Eaux de drainage":
    return  "<CLR green='168' blue='132' >" +[ALTITUDE_FIL_EAU_M] + "</CLR>"
  elif [CONTENU] == "Eaux pluviales":
    return  "<CLR green='77' blue='168' >" + [NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] + "</CLR>"

I need some help finding the right syntax, because i have this error message (image).

If I write code as below, only labels "Eaux pluviales" are displayed

def FindLabel  ([CONTENU], [NOM_OBJET_USUEL], [ALTITUDE_FIL_EAU_M] ):
  if [CONTENU] == "Eaux usées":
    return  "<CLR red='168' >" +[NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] +"</CLR>"
  elif [CONTENU] == "Eaux pluviales":
    return  "<CLR green='77' blue='168' >" +[NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] + "</CLR>"
  elif [CONTENU] == "Eaux mélangées":
    return  "<CLR red='79' green='41' >" + [NOM_OBJET_USUEL] +'\n' + [ALTITUDE_FIL_EAU_M] + "</CLR>"
  elif [CONTENU] == "Eaux de drainage":
    return  "<CLR green='168' blue='132' >" +[ALTITUDE_FIL_EAU_M] + "</CLR>"


texts:
I am trying to format a date field that underwent incorrect importation.  The result was that the date format became coded in a way i am unsure of.

Is there a simple way using the field calculator to batch process date formats?

In the meantime, I will try to create a unique ID in order to join the two tables.

Here is what it currently is:

This is the format I want it in:


texts: to extract data from tif files using QGIS

texts:
I am brand new to GIS. I haven't been able to find the solution online, I think I don't know enough to be searching with the right terminology.

I have a number of tif files which I have successfully uploaded to QGIS as a raster layer. I also have a shapefile which I have uploaded as a vector layer. I am interested in getting the average pixel value for all those pixels of the tif file(s) which correspond to the area delimited by the shapefile.

I seek a solution I can perform in the Python console.


texts: QGIS quits unexpectedly after install "Import Photo" plugin on MX Linux

texts:
I install "import photo" plugin in QGIS 2.14 in MX Linux 18 via menu Plugins >  Manage and Install Plugin. No message on success or not, then QGIS quits unexpectedly. Then I try to open QGIS again. But, QGIS quits again. So, I uninstall QGIS and install QGIS again. But, QGIS still quits unexpectedly.

Then I update and upgrade to MX Linux 19. QGIS still quits unexpectedly. So, I uninstall and install QGIS. QGIS still quits unexpectadly. So, I delete file user.py and user.pyc in /usr/lib/python2.7/dist-packages/qgis/. There is a message:

  Couldn't load qgis.user. Python support will be disabled.

  Traceback (most recent call last):   File "", line 1, in    File
  "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level) ImportError: No module named user

Then I close. After that QGIS runs well. But, I can not install plugin anymore because in Setting tab (Plugins > Manage and Install Plugin > Setting) there is message:

  "The settings on this tab are only for Python Plugins. No Python
  support detected, thus no setting available."

So, any suggestion?


texts: GRASS r.in.lidar problem on Windows 10: 'no such file or directory'

texts:
I'm having problems to load .las files into Grass GIS 7.0.3 for Windows 10. I use the module r.in.lidar but I get the following message:

What can I do? How can I solve the problem?


texts: Extract coordinates from vector layer in PyQGIS

texts:
In QGIS, I have made a polygon, and from that polygon I have made a vector grid via research tools-> vector grid. The grid covers the entire polygon.

Now I want to extract these grid points in the python console so I can work with them. But I can not figure out how to extract these coordinates.

I have done the following:

Layer = qgis.utils.iface.activeLayer()
provider = Layer.dataProvider()
feat = QgsFeature()
allAttrs = provider.attributeIndexes()
provider.select(allAttrs)
geom = feat.geometry()

but from here I do not know what to do. Any help?


texts: QGIS Python matching list values to feature ID

texts:
Working in QGIS 3 x64 with built-in python console and editor.

I have a list of the feature IDs of duplicate features in my dataset. Ultimately, I want each of these features to have a value of 1 for the field error.

When I run the following block of code, the correct number of duplicate features prints, so that's not a problem. It's when I go to update error to 1.

Code so far:

geoms = dict()
dup_features = list()
null_features = set()

features = layer.getFeatures()
index = QgsSpatialIndex()

for current, f in enumerate(features):
    if not f.hasGeometry():
        null_features.add(f.id())
        continue

    geoms[f.id()] = f.geometry()
    index.addFeature(f)

unique_features = dict(geoms)

for feature_id, geometry in geoms.items():
    if feature_id not in unique_features:
        continue

    candidates = index.intersects(geometry.boundingBox())
    candidates.remove(feature_id)

    for candidate_id in candidates:
        if geometry.isGeosEqual(geoms[candidate_id]):
            dup_features.append(feature_id)

print(len(dup_features))

## problem section of code

for feature in dup_features:
    feature['error'] = 1 ## Line 58 referred in error
    layer.updateFeature(feature)

Error generated:

    File "<string>", line 58, in <module>
TypeError: 'int' object does not support item assignment

Alternatives Tried: (none of these throw errors but don't update the attribute value)

for feature_id in dup_features:
    layer.getFeature(feature_id)
    feature['error'] = 1
    layer.updateFeature(feature)

for feature_id in dup_features:
    feature['error'] = 1
    layer.updateFeature(feature)


texts: How to perform Noise Pollution Mapping using ArcGIS for Desktop?

texts:
I've been trying to do a noise pollution map, and I really don't know where to start with the information I have. The end result needs to be like in the link below.

I've seen the question, Creating noise maps using ArcGIS for Desktop or Open Source Desktop GIS? but it doesn't give an accurate metod, and I don't know how to go about it. I've searched for weeks and weeks now about a simple method, and I'm none the wiser.

Data I have:

parameter - speed (roads)
environments (hard ground=lakes and soft ground= grassland)
sound data i.e. dB accounting to number of vehicles

What are the tools I can use in ArcGIS for Desktop to achieve a map like the link above?

I only have access to ArcGIS for Desktop.


texts:
I am trying to spatially thin occurrence data for use in a Species Distribution Model. Because my data represent "presence-only" or "citizen science" data, I am trying to account for the effects of sampling bias by subsampling the points so that no two points fall within X distance of each other. This is the Systematic sampling approach detailed in Fourcade etal. 2014 and many other recent Species Distribution models.
So far, I have tried variations on this theme:
Create spatial buffer around each occurrence point:
var hydrilla = users/me/myoccurrencecsv.csv
var hyd_buffered = hydrilla.map(function(feat){
  return feat.buffer(5000)
})

//Join data so that no two points fall within 5 km

var distFilter =  ee.Filter.intersects({
  leftField: '.geo',
  rightField: '.geo',
  maxError: 100
});

var join = ee.Join.inner()
var hydrilla_corrected = join.apply(hydrilla,hyd_buffered,distFilter)

I am still fairly unfamiliar with the ee.Join process, so I'm sure something I am doing is wrong, but I have tried an ee.Join.inverted() as well as an ee.Join.simple() but nothing seems to give me what I want. I have done an extensive googling of this problem in EE and have not found any satisfactory results.
Edit: Here is another pseudo solution that I propose
var hyd_corrected = hydrilla.map(function(feat){
//Buffer around each feature
  var buffered = feat.buffer(5000)

//Cast each buffer object to a geometry object by referencing the geometry properties of the feature
// and receiving the coordinates
  var geom = ee.Geometry.Polygon(ee.List(buffered.geometry().coordinates()));

//Return Null if the feature intersects the geometry
  return ee.Algorithms.If(feat.intersects(geom),null,feat)
})

This doesn't work because each feature obviously intersects its own buffered geometry, but what is relevant is adjacent features.
Here is the link to my script: https://code.earthengine.google.com/190c7ec503bd9aff08099c7b500ad8a6


texts: Meridian convergence. True north, magnetic north and grid north

texts:
I'm having difficulty wrapping my mind around meridian convergence, the horizontal angular distance between true north and grid north.

Say I want to determine if a feature lies to the east or west of the central meridian (CM) by using true north and grid north designations on a topo map.

I realize grid north is parallel to the CM, while true north points to the poles, and magnetic north is the direction a compass points and constantly changing.

Looking at the figure below, I see grid north points to the east of true north. So does this fact alone mean everything on the grid is to the east of the true north, and therefore the CM? Does this simple fact mean that features in this grid lie to the east of the CM, or am I missing something.

Also, I realize I can just look at the tick marks and and see that this this grid in particular covers -80°22’ to -80°15, and the CM for zone 17R lies at 81°. But had I not known that and wanted to use the north arrows. Please explain.


texts: QGIS select features from attribute table using the "LIKE" operator

texts:
I created an expression with the help of the GUI in the QGIS attribute table of a vector (shapefile) layer of Native American ("Indian") reservations in the USA state of New Mexico.

"NAME"  LIKE 'Picuris'

No row became selected.

The full "NAME" in the desired row is 'Picuris Indian Reservation'  . I don't mind typing the whole name, but in principle, is not the LIKE function supposed to find partial matches? That is, if I want 'State Hwy 68' (which is in the FULLNAME field) (of a different layer), will the LIKE operator not find it if I use "FULLNAME" LIKE 'Hwy 68' ?


texts: Creating new layer based on existing layer Google Earth Engine

texts:
I have many layers consist of many parameters. I want to create a new layer called 'Reforestation Priority Area' based on weighting given to existing layers. Is there any way to do that?
The code is below:
var countries = Indonesia

var Nama_Provinsi = ['Kalimantan Selatan']
var geometry = countries.filter(ee.Filter.inList('NAME_1', Nama_Provinsi));
var area1 = ee.FeatureCollection(geometry);

var image = ee.Image('WCMC/biomass_carbon_density/v1_0/2010');

Map.addLayer(
    image.clip(area1), {
      min: 1,
      max: 180,
      palette: ['d9f0a3', 'addd8e', '78c679', '41ab5d', '238443', '005a32']
    },
    'carbon_tonnes_per_ha');

var dataset = ee.Image('UMD/hansen/global_forest_change_2020_v1_8');
var treeCoverVisParam = {
  bands: ['treecover2000'],
  min: 0,
  max: 100,
  palette: ['black', 'green']
};
Map.addLayer(dataset.clip(area1), treeCoverVisParam, 'tree cover');

var treeLossVisParam = {
  bands: ['lossyear'],
  min: 0,
  max: 18,
  palette: ['yellow', 'red']
};
Map.addLayer(dataset.clip(area1), treeLossVisParam, 'tree loss year')

var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2020-01-01', '2020-12-31'));
var evi = dataset.select('EVI');
var eviVis = {
  min: 0.0,
  max: 8000.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
Map.setCenter(116, -3, 8);

var table_bound = function(image){
  return image.clip(area1);
};
var eviclipped = evi.map(table_bound)

Map.addLayer(eviclipped, eviVis, 'EVI2020');

var dataset = ee.Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019")
.select('discrete_classification');

Map.setCenter(116, -3, 8);

Map.addLayer(dataset.clip(area1), {}, "Land Cover");

//displaying the data in interactive map
Map.addLayer(image.clip(area1), {min:0, max: 300,
palette: ['90EE90','FFFF00','FF0000']}, 'Raw SRTM'); //Light green, yellow, red

//display hillshading and slope
var hillshade = ee.Terrain.hillshade(image);
Map.addLayer(hillshade.clip(area1), {min:150, max:255,}, 'Hillshade', false);
print(hillshade);

var slope = ee.Terrain.slope(image);
Map.addLayer(slope.clip(area1), {min:0, max:20, pallete: ['FFFFFF']},'Slope', false);
print(slope);


texts: Animation does not work GeoJSON date fields in CARTO VL

texts:
In CARTO VL linear function expects a date or numeric field as its first argument. But when the data source is a GeoJSON, date fields are stored as categorical/strings and it returns an error like this:

CartoValidationError: [Property with an incorrect type]: linear(): invalid first parameter 'input' expected type was one of number,date, actual type was 'category'

You can test it in this visualizer replacing the current CARTO VL viz style with something like:

 filter: animation(
             linear(
                    $date,
                    time('2019-01-01T00:00:00Z'),
                    time('2019-01-31T00:00:00Z')
                ),
                20,
                fade(0.5, 0.5))

Is there any workaround to convert this (string/categorical) date field to a proper date/timestamp field in CARTO VL?


texts: Why is Shapely and contains_point not detecting point being on land?

texts:
I'm using the world shapefile that can be found at https://trac.openstreetmap.org/browser/applications/editors/merkaartor-branches/merkaartor-0.14-fixes/share/world_shp/world_adm0.shp?rev=17052.

I've plotted the point (150.625, -26.7396) (longitude, latitude) using QGIS and clearly it is on land.

Here's my code:

import pysal as ps # asShape, contains_point;
import shapely.geometry # asShape, buffer

WORLD_FILE = 'world_adm0.shp'
BUFF = 0.005

polygons = ps.open(WORLD_FILE)

found = False

for polygon in polygons:
    sh_poly = shapely.geometry.asShape(polygon)
    sh_buff_poly = sh_poly.buffer(BUFF)
    buff_poly = ps.cg.asShape(sh_buff_poly)

    if buff_poly.contains_point((150.625, -26.7396)):
        found = True
        break

if found:
    print("on land")
else:
    print("in sea")

Why do I get "in sea"?


texts: Assigning field value based on location without join - Two different feature layers

texts:
I'm interested in updating a field in N-geodatabase feature classes based on their location to another feature class in the geodatabase, but without having to perform join.

I am new to Python, and have tried to use the similar answer to my question in : Assigning field value based on location without join

I have a feature layer "zone - M" and a feature layer "Zone - Z", these features have a field called "name_M" and "name_Z".
I will draw several others layers "on" these zones. These others layers have all a field "name_zone_M" and "name_zone_Z". I would like that if the layer is contained completely inside a zone, the value "name" of each zone get assigned to the field "name_zone_M" and "name_zone_Z".
These feature class are all polygons.

I have tried the following code to try to get the first one right, and then to try to loop inside the others layers, but I did not succeed.

Any toughts?

import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
# Create FeatureLayers
arcpy.MakeFeatureLayer_management(r"Zone-M")
arcpy.MakeFeatureLayer_management(r"DRAW1")

# Create a search cursor for the states
rows = arcpy.SearchCursor(r"Zone-M")
for row in rows:
    # What you'll do is select each neighborhood one at a time, and then select all the parcels in that neighborhood and calculate the SITE_CODE field
    # NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
    arcpy.SelectLayerByAttribute_management(r"Zone-M", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
    arcpy.SelectLayerByLocation_management(r"DRAW1", "COMPLETELY_CONTAINS", r"Zone-M", "", "NEW_SELECTION")
    arcpy.CalculateField_management(r"DRAW1", "name_zone_M", "'{0}'".format(str(row.getValue("name_M"))), "PYTHON_9.3", "")
    print "Finished processing " + str(row.getValue("name_M"))


texts: Obtaining NDVI statics per region and year using Google Earth Engine

texts:
How do I obtain a table which collects the mean NDVI per year of the region like this?
2005: 0.37
2006: 0.31
2007: 0.33
2008: 0.40
2009: 0,35

I have this code; I need to finish it and I don't know how.

//seleccion del pais
function filtro_pais (pais_eleg) {
  var filtro = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
    .filter(ee.Filter.or(
        ee.Filter.eq('country_na', pais_eleg)));

    return filtro;
}
var pais_elegido = "Spain";
var fc = filtro_pais (pais_elegido);

//invocacion a la base de datos de Modis
var imagen1 = ee.ImageCollection('MODIS/MCD43A4_006_NDVI')
                .filter(ee.Filter.date('2015-01-01', '2016-01-01'))
                .select('NDVI');

 var mean1 = imagen1.mean();

  var clip1 = sum1.clipToCollection(fc);

  var sum = clip1.reduce('mean')


texts:
I have a simple application that uses OpenLayers 2.13.1 to draw GML data pasted in a text box. There is already some default GML preloaded in textarea and it works fine. When you press "Add feature" button, it draws polygons on the map.

GitHub link:
https://github.com/mleontenko/OpenLayers2-DrawGML/blob/develop/index.html

Now, I have some other GML data that i get with WMS GetFeatureInfo request. It looks like this: https://raw.githubusercontent.com/mleontenko/gml/master/croatia.gml

But this new GML doesn't draw visible polygons on map. Does anyone know why? And, what do i need to do to draw this GML correctly?


texts:
Always when working with raster layers errors happen. Can´t warp (reproject) layer, or extract contour, for example.
I have unistall Qgis and install back several times in different ways, but didn´t work.
I work with MacBook. This is an example of one of the errors:

  The following layers were not correctly
  generated./var/folders/0z/1728_nfx30j9x0vgs79sb9zr0000gn/T/processing_b2c1d0382bbe42b9936283783605c8a7/840b94611d8347b2936535f235e21e98/OUTPUT.tif/var/folders/0z/1728_nfx30j9x0vgs79sb9zr0000gn/T/processing_b2c1d0382bbe42b9936283783605c8a7/a410e5fa0c8444c88370c60c91b55fff/OUTPUT.tifYou
  can check the 'Log Messages Panel' in QGIS main window to find more
  information about the execution of the algorithm.


texts: How to create raster layer directly from numpy array in QGIS, without saving to file?

texts: I have managed to create a numpy array from a raster file using commands from the gdal library, and now I would like to load that raster as a layer without saving it first. It seems like there should be a simple way to do this. Is there?


texts:
I'm somewhat familiar with ArcGIS Desktop, but have no experience in Python. I have about 100 shapefiles which have an attribute table containing a date, but I need a column which contains a corresponding season.

To complicate things a little more, I do not want to define seasons the traditional way. Seasons are defined as follows:
a.  Summer = June, July, August, September
b.  Fall = October, November
c.  Winter = December, January, February, March
d.  Spring = April, May

So far, I've been pasting and copying each attribute table into Excel, and using a function that looks like this: =CHOOSE(MONTH(J2),"Winter","Winter","Winter","Spring","Spring","Summer","Summer","Summer","Summer","Fall","Fall","Winter")

I'm sure there is a faster way to do things.

I recently added a 'MONTH' column using this script: Left( [GMT_DATE] , Len( [GMT_DATE] ) -8 )

I figured out how to generate a season column from a month column, but I am still unable to generate a season column directly from a date column. When I try I get "Winter" for everything.


texts:
So initially I had multiple perpendicular lines along a coastline.

Next I clipped these lines based on polygons along the coastline,such as the follwing.

I am trying to use this resulting shapefile to get the pixel values from a raster based on these clipped lines. The problem is that I am not being able to iterate over these as I would before I clipped them.

This is the script I'm trying to use

from osgeo import ogr, gdal
import rasterio
import numpy as np
from shapely.geometry import MultiLineString, LineString, Point
from shapely import wkt

raster = rasterio.open(r'C:\Users\Cassiano\DATA_ANALYSIS\MAPBIOMAS\1985\mapbiomas-riograndedosul-1985.tif')
band = raster.read(1)
ds = ogr.Open(r'C:\Users\Cassiano\DATA_ANALYSIS\MAPBIOMAS\LINHAS\RS.shp')
layer = ds.GetLayer(0)

#Setting lists to store data
x = []
y = []
z = []
line_id = []
# distance of the topographic profile
distance = []
dist = 0.03 / 111.139

This is the snippet that really matters

for line in layer:
        line_points = []
        cur_dist = dist
        line_geom = line.geometry().ExportToWkt()
        shapely_line = LineString(wkt.loads(line_geom))
        line_len = shapely_line.length
        while cur_dist < line_len:
            point = shapely_line.interpolate(cur_dist)
            xp, yp = point.x, point.y
            line_id.append(line.GetFID())
            x.append(xp)
            y.append(yp)
            #get pixel value
            row, col = raster.index(xp, yp)
            z.append(band[row, col])
            distance.append(cur_dist)
            cur_dist += dist
profile = zip(line_id, x, y, z, distance)

This is the error message

NotImplementedError                       Traceback (most recent call last)
<ipython-input-19-2f8fa74d07bc> in <module>()
      3         cur_dist = dist
      4         line_geom = line.geometry().ExportToWkt()
----> 5         shapely_line = LineString(wkt.loads(line_geom))
      6         line_len = shapely_line.length
      7         while cur_dist < line_len:

~\Anaconda3\lib\site-packages\shapely\geometry\linestring.py in __init__(self, coordinates)
     46         BaseGeometry.__init__(self)
     47         if coordinates is not None:
---> 48             self._set_coords(coordinates)
     49
     50     @property

~\Anaconda3\lib\site-packages\shapely\geometry\linestring.py in _set_coords(self, coordinates)
     95     def _set_coords(self, coordinates):
     96         self.empty()
---> 97         ret = geos_linestring_from_py(coordinates)
     98         if ret is not None:
     99             self._geom, self._ndim = ret

~\Anaconda3\lib\site-packages\shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linestring_from_py()

~\Anaconda3\lib\site-packages\shapely\geometry\base.py in __array_interface__(self)
    792     def __array_interface__(self):
    793         """Provide the Numpy array protocol."""
--> 794         raise NotImplementedError("Multi-part geometries do not themselves "
    795                                   "provide the array interface")
    796

NotImplementedError: Multi-part geometries do not themselves provide the array interface

In order to get the lines clipped I used the clip method from QGIS with the lines and the polygons.
This script has worked for other shapefiles, I started getting this error message only after clipping with QGIS.

I decided to use QGIS because I was having some trouble clipping the lines by writing my own code, though I don't understand what could've compromised the geometries since whenever I try the argument

line.geometry().GetGeometryName()

The result is always a 'LineString' or a 'MultiLineString', therefore iterating through these features shouldn't be a problem... I guess.


texts: Buffer to linestring in shapely produces extra line


texts:
I have given this code for intersection:
function intersectPolyByPolyFC(poly,fcPoly){
        var fgp=[];//variable for feature group
        var bbPoly=turf.bboxPolygon(turf.bbox(poly));//polygon representation of the bounding box of a reference polygon
        for(var i=0;i<fcPoly.features.length;i++){//loop through the feature collection of polygons
            var bb=turf.bboxPolygon(turf.bbox(fcPoly.features[i]));//variable bb represents each polygon for feature collection
            if(turf.intersect(bbPoly,bb)){//condition to check if the two bounding boxes intersect
                var int= turf.intersect(poly,fcPoly.features[i]);
                if(int){
                    int.properties=fcPoly.features[i].properties;
                    fgp.push(int);
                }
            }
        }
        return turf.featureCollection(fgp);
    }

and the function definition I have given as
var bufline=turf.buffer(line,0.05,{units:"kilometers"});
var intBUOWL= intersectPolyByPolyFC(bufline, lyrBUOWL.toGeoJSON());
L.geoJSON(intBUOWL,{style:{color:'red',weight:5}}).addTo(mymap);

When I am drawing a line over a polygon for intersection I am getting this error
Uncaught TypeError: Cannot read property 'map' of undefined
    at Object.t.collect (turf.min.js:1)
    at summarizePointsByLine (buffer.html:440)
    at NewClass.<anonymous> (buffer.html:155)
    at NewClass.fire (Events.js:190)
    at NewClass._fireCreatedEvent (leaflet.draw.js:8)
    at NewClass._fireCreatedEvent (leaflet.draw.js:8)
    at NewClass.completeShape (leaflet.draw.js:8)
    at HTMLAnchorElement.handler (DomEvent.js:92)

How can I resolve this error?


texts:
All,
I have 2 layers as shown below

What I want to do is go through each feature (in pink) and do a spatial filter against the green layer.
My current code is

        inMemFeatCur = inMemFeatClass.Search(null, true);
        inMemFeature = inMemFeatCur.NextFeature();
        while (inMemFeature != null)
        {
            //set the spatial filter
            ISpatialFilter spaFil = new SpatialFilterClass();
            spaFil.Geometry = inMemFeature.Shape;
            spaFil.GeometryField = validationLayer.ShapeFieldName;
            spaFil.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
            spaFil.SpatialRelDescription = "";

            //query the results from the spatial filter

            IFeatureCursor geomCursor = inMemFeatClass.Search(spaFil, true);

The problem is the when I go through the geomCursor, it will not find the polygon that is completely outside.

The check i want to do is each polygon (pink) completely within (common boundaries allowed) the green layer. If it is not, a warning will to appear.

I know that IRelationalOperator returns a bool, but i have to go through each green feature which I don't want to do, as I'll be looping many times for each pink feature.

Has anyone got any ideas on how to do this.


texts: How to mask only select layers (i.e., exclude some layers from the mask)

texts:
is there a way to create a masking layer in QGIS (2.14) that will mask only certain layers, so that some layers (such as baselayer imagery or other) will be unaffected?

I'm using geometry generator as a symbol layer type to create dynamic, virtual buffers and want them to be masked when outside of a particular feature such as a property boundary

I have tried the mask plugin and also tried the inverted polygon symbol on a layer but cannot achieve desired effect as I can't find an option to exclude or include layers.


texts:
I have layer of buildings (polygons), and in attribute table contains number of buildings, floor area of buildings etc. I would like to know how many buildings (#) is per Ha or km2, and building area per Ha or km2.

How can I calculate that in QGIS3?


texts: How to run the squarring algorithm from opencarto on a shapefile with many polygons?

texts:
I want to square building surfaces which are almost rectangular. I didn't find any tool in an OSgis. But there is a java library opencarto which includes such an algorithm. I'm not a professional java developer, but as i can see in the code below there are two getter methods which wait for a polygon input.
My Question is how can I select with a Filechooser in java the shapefile and make every single polygon of this shapefile running through this algorithm?
Or exists there almost a tool which can handle this?

public class Squarring {

public static Polygon get(Polygon poly) {
    return get(poly, 20, -1, 0.01);
}

public static Polygon get(Polygon poly, double angleToleranceDeg, int pointActivationLimitNb, double res) {
    Polygon poly2 = (Polygon) poly.clone();
    HashSet<GPoint> ps = new HashSet<GPoint>();
    HashSet<GSegment> segs = new HashSet<GSegment>();
    Decomposers.decomposeLimit(poly2, res, ps, segs, false, null);

    double o=new Orientation(poly2).getSidesOrientation();

    for (GSegment s : segs)..........


texts: Which version of GQIS should I download Windows 7 64 bit PC?

texts: I have downloaded GQIS 2.4 on my PC and have gone through the installation but each time I try to open download, it fails and then prompts to try to reinstall it. I am unsure if I am using the correct version of the download I need and how to correctly install it to my PC.


texts:
I'm trying to set up a moving map on my tablet (asus, Android 4).
I've installed QGIS (latest version), which works. I've got a blue tooth external GPS, which I've paired to the tablet. However, I cannot link the GPS to QGIS.

I followed the steps on QGIS website "go to View-> Panels-> " but no other page comes up.

Any suggestions?


texts: Cutting lots of polylines using single horizontal line in ArcGIS for Desktop?

texts:
I have 2000 polylines or more, each polyline have about 2000-3000 vertices. I want cut all part of polylines that are below a single horizontal line. How can I do that in Arcgis?

I try to use intersect tool: polylines layer with the layer horizontal line (with output "point"). In this way I obtain points intersection, I can use it to split the polylines and finally I can delete all polylines below the horizontal line by simply selection the area. The problem is that the first step, intersecate polylines with the single horizontal line, take an infinte process. My idea isn't good, is there another fast solution?


texts: Trying to fetch feature data by intersecting Vector Polygon and WMS feature layer

texts:
I am trying to intersect a user defined polygon with a Point WMS Feature Layer. I am using the following code

action = new GeoExt.Action({
        text: "draw poly",
        control: new OpenLayers.Control.DrawFeature(
            features, OpenLayers.Handler.Polygon, {
            eventListeners: {
                "featureadded": function(evt) {
                    var pfilter = new OpenLayers.Filter.Spatial({
                             type: OpenLayers.Filter.Spatial.INTERSECTS,
                             value: evt.feature.geometry
                         });

                }
            }
        }),
        map: map,
        // button options
        toggleGroup: "draw",
        allowDepress: false,
        tooltip: "draw polygon",
        // check item options
        group: "draw"
    });
    actions["draw_poly"] = action;
    toolbarItems.push(action);

Now I do have the polygon geometry,

POLYGON((-12440279.225738 3759478.9686547,-12376683.618213 3731350.1422497,-12429272.293666 3708113.2856542,-12440279.225738 3759478.9686547)) { id=

"OpenLayers_Geometry_Polygon_384"

, components=[1], bounds=-12440279.225738,3708113.2856542,-12376683.618213,3759478.9686547, more...}

But I am not sure how do I proceed from here? Pls help.


texts: Adding a custom color scale for a sp.points layer to levelplot() legend

texts:
I'm wondering how I can add a customized color palette that I used to display a sp.points layer in a call to levelplot (lattice plotting) in the rasterVis package in R.

I used a function that Dave X created here in his answer to a question on how to map numeric values to colors. Now, I am using it in side-by-side maps with different raster backgrounds, as the color scale that a points layer is mapped to. I want this custom color ramp to appear in a legend but I don't need it to show up twice - it's the same on both panels of the plot.

A quick summary of the maps I'm making: A two-paneled plot of

elevation, and
aspect (calculated from the elevation layer),

for a district in rural Tanzania. Plotted over these rasters on each panel is a polygons layer denoting the district's adm boundaries, and a points layer representing grain mills in the district that I visited on a prior field trip. The color represents the mean Aflatoxin level observed in samples of milled maize taken from each mill. The code for one of the panels of the plot (to give you an idea) is below:

(See Dave X's post linked above for the code for map2color() )

    AF <- MillGeo@data$AF_mean
    AF_colors <- colorRampPalette(c("light yellow", "dark red"))(length(AF))
    pal <- map2color(AF, AF_colors)
    p1 <- levelplot(elev_raster, margin=FALSE,
            at=seq(600,2200, by=50), col.regions=terrain.colors(160),
            xlab=NULL, ylab=NULL) + layer(sp.polygons(Kongwa)) +
            sp.points(MillGeo, pch=20, cex=1.5, col=pal)

Imagine that there are two such plots next to one another and they both use the same color palette for the points representing Aflatoxin levels. I want to customize a legend that I can represent the scale according to the values of Aflatoxin used to create it, ideally just showing up once below the two panels. Is there a way to do this?

Right now, my code for showing the plots is:

    grid.arrange(p1,p2, ncol=2)


texts: Is there a clever way to combine two or more geoJSON files? I would like to combine the two geoJSON files here: https://gist.github.com/szeconku/3c0d3021178d8598f04dc5273864bd81


texts:
I have a shapefile 1 and a shapefile 2. Shapefile 1 contains ~2k Polygons while shapefile 2 contains ~4k polygons.
I want to know, which Polygons from shapefile 1 intersect polygons in shapefile 2. Further want to write the information from the intersecting polygons into the attributetable from shapefile 2, into the corresponding rows of the polygon. I want to preserve the shape/outline of shapefile 2!
One polygon from shapefile 2 can be intersected by many polygons from shapefile 1. That means my resulting shapefile will have way more than ~4k Polygons. Whatever I try though (Intersect/Union) returns me the information I need BUT with the outline of shapefile 1. I tried a Spatial Join but that did not help neither. It did not write multiple information into one polygon in shapefile 1. The resulting shapefile from the spatial join had as many rows as shapefile 2.

Is there a way to do this in QGIS/GRASS/GDAL?


texts: stack rasters with different extent and resolution in a loop with R package raster

texts:
I would like to stack rasters with the R package "raster". Layers are in ascii format. I want to crop them to have the same extent and project them to have the same resolution. My code gives raster layers with empty values.

#The layers must have a similar extent and projection. We will first crop them all to the same extent.
#Read files found in the working directory containing *.asc
raster_data <- list.files(pattern='*.asc', full.names=TRUE)
#get desired extent (smallest raster extent)
e <- extent(-12, 112, 23, 72)
#get template layer with desired extent and resolution. Transform ascii template layer to raster format
d <- raster("WC_BIO1.asc", ext = e, crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0", ncol=15120, nrow=6000)
w <- writeRaster(d, "WC_BIO1_template.grd", format = "raster")
#load ascii layers, transform to raster format with extent and resolution obtained from template layer
for (i in 1:length(raster_data)){
  r <- raster(raster_data[i], crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  rc <- crop(r, e)
  rp <- projectRaster(from = rc, to = w,
                      crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
                      filename = file.path ("./crop", raster_data[i]),
                      method = "bilinear",
                      format = "raster",
                      overwrite = TRUE)
  rw <- writeRaster(rp,
                    filename = file.path ("./crop", raster_data[i]),
                    crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
                    overwrite = TRUE,
                    format = "raster",
                    na.value=-9999)
  }


texts: Calculate mean value from list features in Google Earth Engine

texts:
I created a FeatureCollection containing 63 features with a different date, from one filter per ID.

In each of these 63 features, I created a new property named "matches" that contains a new list of features, which stores a value (called LANDSAT) for each of the 50 samples I disperse in my region.

Now, I would like to redeem the LANDSAT value of each list item and take the average, how would I do this?

 var collection = Collection .select(["NDVI_to_ajust"], ["LANDSAT"])

      var triplets = collection.map(
        function(image) {
        return image.reduceRegions({
          collection: samplesVegPrimaCinquenta,
          reducer: ee.Reducer.mean().setOutputs(image.bandNames()),
          scale: 30,
        }).map(
          function(feature) {

          return feature.set({
            'imageID': image.id(),
            'date': ee.Date(image.get('system:time_start')).format('Y-MM-dd')
          })
        });
      }).flatten();
      print(triplets,'triplets')

    var coord = function(table,rowId, colId, rowProperty, colProperty){
      var rows = table.distinct(rowProperty);
      var joined = ee.Join.saveAll('matches').apply({
        primary: rows,
        secondary: table,
        condition: ee.Filter.equals({
          leftField: rowProperty,
          rightField: rowProperty
        })
      });
    print(joined,'join by same ID/DATE')

Here's the link code


texts: Arcgis Javascript API WMTS Layer don't use the configured proxy

texts:
I am new about Arcgis javascript API. I'm using WMTSLayer and i have configured to use a proxy.

esriConfig.defaults.io.proxyUrl = <Proxy_url>

When the API does the GetCapabilities request, the request goes through the proxy. But tiles requests does not use the proxy.

This is the expected result? What can i do to force the tiles request to use the proxy configuration?


texts: How to Clear Selected Features in Arcmap using Python

texts:
I would like to Clear Any Selected Features using python since it causes an error if I run my tool and the user has a recorded selected.

I need some thing here that runs ("CLEAR_SELECTION") first

rows = arcpy.SearchCursor(Table, " \"Name\" LIKE 'Albert%' ", "", "Name;")

Thoughts?


texts:
Anyone know what I can use to convert a geoid surface stored as , BYN into .GTX for use with PROJ4?

The CCG2013A vertical datum and geoid model from NRCan is only available, as far as I can find, in BYN format.


texts: Rendering point symbols on top of diagrams?

texts:
I've got some data where I have categorized symbology for point data with different colors showing different categories. I would also like to display some data using the pie charts function. I would like to have the different colored point on top of the pie charts to get a symbol that looks something like the mock-up below.

I can't find a way to get the points on top of the charts. Any suggestion?


texts: ST_ClosestPoint in 2D but retaining Z value

texts:
I have a query in PostGIS where I have a two geometries, a MultiLineStringZ and a Point.

I want to find the 3D point on the MultiLineStringZ closest to the Point by it's flat 2D plane distance, so essentially ST_ClosestPoint(myMultiLineStringZ, myPoint). However I need to retain the Z value on the returned point. ST_ClosestPoint drops the Z, and ST_3DClosestPoint does a 3D distance calculation.

How do I do this?


texts: Cartodb image header doesn't show up after publishing the map

texts:
I've linked image header to my map, I normally can view it while editing it. Then, when I publish it via link or embed code, the image doesn't appear, only a gray color.

I customized HTML code only by changing the height. This is the code, is there something wrong?

<div class="cartodb-popup header with-image v2" data-cover="true">

x

        Non-valid picture URL

      նկարագրություն
      {{description}}


texts: Points Disappearing When Applying Categorised Symbology in QGIS

texts:
I am new to QGIS (and GIS in general).
I have imported a CSV file containing approx 1500 points with lat/long co-ordinates. I applied categorised symbology based on a field in the attribute table, and that worked with no issues, so I now have points with various colours based on the values in a field:

I am now trying to amend some of the data in the attribute table so that I can move some points from one category into another category.
As this layer is based on a .csv file, I am unable to edit the attribute table. I read some posts where people had recommended saving a copy of the layer as a Shapefile or a Geopackage, to allow editing of the attribute table.
So I have done this, but this removes the original categorised symbology and all of the points are back to one colour in the new layer:

I then try to re-apply the categorised symbology to this new layer as follows:

...but as soon as I click 'Apply' or 'OK', the points disappear from the map:

I am not sure what the issue is or what I am doing wrong, does anybody have any ideas?


texts: Cannot write a subset of features back to the source Geopackage or File Geodatabase

texts:
At QGIS 2.18.16, I have a Geopackage point layer with a selection subset.  I need to write the selected records to the same Geopackage as a new layer, with a new name.

When I attempt to Save As... the new layer, QGIS responds with this error:

It appears that QGIS thinks that I'm trying to overwrite the original layer, when I'm actually trying to create a new layer that is a subset of the original.

Interestingly, the same error occurs if I attempt the same with a File Geodatabase; that is, I try to write a subset of records as a new Feature Class to the original FGDB.

A Google search turned up the same error, but in that case the user was trying to overwrite the source layer.

I expect that QGIS should be able to write a subset of a layer's records to the source Geopackage or File Geodatabase under a new name.  How can I accomplish this without the error?

As a workaround, I can write the subset to a different Geopackage (or FGDB), and then copy that to the original Geopackage (or FGDB).


texts: qgis2web display only from a certain scale

texts:
I'm using QGIS 3.4 and the latest qgis2web plugin.
I have a problem to display features at on some scales.
Openstreetmap is displayed good but the polygons of a layer are only shown when I zoom in to scale 1:10.000, at smaller scales the features disappear.
The number of features in the layer are only 373.
Can somebody give me a hint where to look for the cause?


texts: How to show latitude and longitude in selected Geometry in Google Earth Engine

texts:
I want to display the selected boundary with Latitude and Longitude. Please guide me on how to I can do this.
here is my code:
/////////////////////////////////////////////////////Bergamo
var BergamoGeometry = ee.Geometry.Polygon([
 [9.66775774583506,45.68010702771403],
[9.69797014817881,45.68346506712161],
[9.69144701585459,45.708643938278705],
[9.660891290756934,45.70217045421897],
[9.66775774583506,45.68010702771403]
]);

// Create a FeatureCollection from the list and print it.
var BergamoGeometry = ee.FeatureCollection(BergamoGeometry);
print(BergamoGeometry);
Map.centerObject(BergamoGeometry);
Map.addLayer(BergamoGeometry);
//////////////////////////////////////

var dataset = ee.ImageCollection("NASA_USDA/HSL/SMAP_soil_moisture")
                 .filterBounds(BergamoGeometry)
                  .filter(ee.Filter.date('2020-03-30', '2020-04-29'));

var soilMoisture = dataset.select('ssm');
///////////////////////////////
Export.image.toDrive({
  image: soilMoisture,
  description: 'SSM_Bergamo',
  scale: 30,
   fileFormat: 'GeoTiFF',
  folder: 'GEE',
  region: BergamoGeometry
});
var soilMoistureVis = {
  min: 0.0,
  max: 28.0,
  palette: ['0300ff', '418504', 'efff07', 'efff07', 'ff0303'],
};
Map.centerObject(BergamoGeometry, 9);

// Create and print the chart.
print(ui.Chart.image.series(soilMoisture, BergamoGeometry, ee.Reducer.mean(), 30));


texts: Clipping large shapefile in ArcMap and .dbf size remains same?

texts:
I'm working with a large point dataset (700,000+) features which is just over 2gb in size. I wanted to split the dataset to reduce file size. However when I do a clip in ArcMap 10.5 on the shapefile, the .dbf remains the same size. I've tried searching for previous questions and read through several documents to no avail. Any thoughts?

Let me provide further context based on comments/suggestions,
I've tried clipping down to about a 1/3 of the points and only a very insignificant amount of file size for the .dbf was shaved off.

Keep in mind that the file is comprised of about 700,000+ features/rows. To try to select 1/3 of the feature by selecting rows within the attribute table would be a very cumbersome process compared to select tools or other methods.

I've also tried using the selected features method, than created layer from selected features, then export to shapefile. No change in file size.


texts:
I am using ArcMap 10.5.1, and I am having trouble getting my labels to draw somewhere besides the center of the map in layout view. I am using Maplex. The problem is that I have a polygon from a different layer right in the center of the map, which is where I want it, and ArcMap insists on placing township/range labels right on top of it. The only way to keep the labels off of the polygon is the shift the view so that the polygon is no longer in the center. I would rather shift the labels than shift the view.

In this map, "Current Inventory" is the layer I need unobscured by the labels for the township/range layer. I can't seem to tweak the position properties to shift the labels away from the center.


texts: Using Google API 3 with Openlayers Brings POPUP

texts:
Possible Duplicate:
  Google Maps layer copyright popup every time map updated / user input

I recently upgraded from google map api v2 to v3. It was working quite fine previously but now it does work but each time I zoom,pan I get a popup. I had look at many thred and https://getsatisfaction.com/opengeo/topics/how_to_get_rid_of_goggle_layers_map_data_windows one especially that explains that openlayes 2.11 has no problem (2.10 has problem and 2.11 being the latest I assume this). I am using openlayers 2.11 with geoext. How can I solve this problem?


texts: Getting intersection of multiple polygons efficiently in Python

texts:
I would like to get the intersection of multiple polygons. Using Python's shapely package, I can find the intersection of two polygons using the intersection function. Is there a similar efficient function for obtaining the intersection of multiple polygons?

Here is a code snippet to understand what I mean:

from shapely.geometry import Point

coord1 = ( 0,0 )
point1 = Point(coord1)
circle1 = point1.buffer(1)

coord2 = ( 1,1 )
point2 = Point(coord2)
circle2 = point2.buffer(1)

coord3 = ( 1,0 )
point3 = Point(coord3)
circle3 = point3.buffer(1)

An intersection of two circles can be found by circle1.intersection(circle2). I can find the intersection of all three circles by circle1.intersection(circle2).intersection(circle3). However, this approach is not salable to a large number of polygons as it requires increasingly more code. I would like a function that takes an arbitrary number of polygons and returns their intersection.


texts: Embeding GeoNetwork's client in my application

texts: Does anybody knows how to embed my GeoNetwork client in my application so that user will be able to publish their metadata while they are still on the application.


texts: Issue on Point Graphic Order in ArcGIS API for JavaScript

texts:
I am trying to add a point to map like this:

var point = new esri.geometry.Point(19.82, 41.33);
point = esri.geometry.geographicToWebMercator(point);
var symbol = new esri.symbol.PictureMarkerSymbol("https://cdn1.iconfinder.com/data/icons/Map-Markers-Icons-Demo-PNG/128/Map-Marker-Marker-Outside-Chartreuse.png", 32, 32);
var graphic = new esri.Graphic(point, symbol);
var layer = new esri.layers.GraphicsLayer();
layer.add(graphic);
map.addLayer(layer);
map.reorderLayer(layer,1);

however the marker is showing on the back of borders layers! I tried to keep the point at first level using:

map.reorderLayer(layer,1);

but it didn't work! Can you please let me know what I am doing wrong?


texts: Getting field from getGetFeatureInfo in OpenLayers?

texts:
I'm making a getGetFeatureInfoRequest() as following:

this.map.on('singleclick',(evt) => {
  var viewResolution = /** @type {number} */ (this.map.getView().getResolution());
  var url = this.layer.getSource().getGetFeatureInfoUrl(
    evt.coordinate, viewResolution, 'EPSG:25830',
    {'INFO_FORMAT': 'text/html'});
  console.log(url);
  if(url){
    var parser = new GeoJSON();
    $.ajax({
      url: url,
      dataType: 'jsonp',
      jsonpCallback: 'parseResponse'
    }).then(function(response) {
      var result = parser.readFeatures(response);
      console.log("RESP:"+result);
    });
  }

});

When the request is made, I receive an html response, which has a <table> and a column with the desired value I need to get.

Is there any way to force the server to return data in JSON?

If not, how can I get that value?


texts: Issue with v.clean (GRASS command) in QGIS

texts:
I've trying to clean topology in certain shapefiles.

When I use v.clean command from GRASS 7 in QGIS, the results is a shapefile with a different CRS from original shapefile. Also, when I try to ejecute another commands over v.clean output (as buffer, i.e), results are vectors with no features.

For example:

I'm working in ETRS89 / UTM 30N by default and without "on the fly" reprojection (1).

Result of v.clean (green line) is accord with original vector (red lines)(2).

But when I check CRS in properties (3), it doesn't have one. Also, v.clean doesn't save a .qpj file with cleaned shapefile, but save a .prj file.

Checking the arguments of .prj file, it's completely different from original shapefile, project CRS and CRS configuration for new layers (original in top, result in bottom, 4).

Setting projection from properties box doesn't work, I can't ejecute basic geoprocesses over the "cleaned" layer.

I've tried in OS X and Windows QGIS 2.18 version without positive results. Is this a bug of last version or I'm doing something wrong?


texts: Using output Feature Dataset of Import XML Workspace process as parameter for continued processing in same model?

texts:
I am setting up a model for some automated SDE creation.  As of now, I manually import XML Workspace schema to a new SDE database, then run a model where I manually input the created Feature Dataset from the XML import.  I would like to incorporate the Import XML Workspace tool into my model.

How do I "grab" that created Feature Dataset and use it as a variable for the rest of my model processes?


texts: Calculating datetime difference using ArcPy

texts:
How would I calculate the time difference between two fields?
For Example: 3/22/2018 7:08:00 AM - 3/22/2018 7:07:30 AM


texts: JTS ConformingDelaunayTriangulationBuilder, what is the right way to set tolerance?

texts:
I have one 3D polygon and one 2D polygon. I know that 2D polygon is inside 3D polygon (when projected on the plane) and I want to find Z-coordinate for 2D polygon points. The idea is to apply triangulation and get Z coordinate value from 3D polygon.

Here is my code:

def getZAtInteraction(threeDPolygon: Geometry, twoDPolygon: Geometry): List[Long] = {

    val builder = new ConformingDelaunayTriangulationBuilder()
    //builder.setTolerance(0.0000001D) ???
    builder.setSites(threeDPolygon)

    twoDPolygon.getCoordinates.toList.map { coordinate =>
      builder.getSubdivision.locate(coordinate).dest().getZ.toLong
    }

Sometimes I get exception

  com.vividsolutions.jts.triangulate.quadedge.LocateFailureException:
  Locate failed to converge (at edge: LINESTRING( XXX1.97904355
  YYY1.35001210000007, XXX2.9790435 YYY2.3499198)).  Possible causes
  include invalid Subdivision topology or very close sites

Google suggests to set tolerance ConformingDelaunayTriangulationBuilder

Is it good suggestion?
What is the right way to estimate tolerance?


texts: Meaning of statement in GRASS r.horizon manual

texts:
I have a question about this statement in the r.horizon manual:

"If your projection is based on distance (easting and northing), these
too must be in meters. The buffer parameters must be in the same units
as the raster coordinates."

I have a DEM TIF based on a Web Mercator projection for an area of interest at latitude -37.841477522046524 (scale factor ~1.26). The pixels in the TIF have a real world resolution of 0.15m.
Do I need to preprocess the TIF in anyway to satisfy the "these too must be in metres" requirement? Do I need to rescale to a pixel size of 1m? Do I need to obtain transform the TIF into a projection where the scale factor is 1?


texts: Looks like cesium-terrain-builder still have this issue and can't be used for large projects, so I wonder if there any alternatives to cesium-terrain-builder that can produce DEM tiles, can work with large projects and compatible with cesium?


texts: Tutorial on migrating Add-ins from ArcGIS 10.x (ArcObjects) to ArcGIS Pro for Microsoft VB.NET?

texts:
I can't seem to find any comprehensive tutorials on migrating (recreating comparable functionality) of 10.x (ArcObjects) VB.Net Add-ins to ArcGIS PRO SDK Add-ins.

When migrating from Avenue, ESRI offered courses, books and tutorials, etc.
When Add-Ins were introduced (around the time VBA was dropped) in version 10, ESRI had a complete A-Z online tutorial on migrating from VBA/VB6 to VB.Net and C# including snippets in both VB and C#.
Both versions 8 - 9 and 10.X used similar ArcObjects Component Object Model so migration path was straight forward in terms of writing the actual to customize ArcGIS using ArcObjects.

I understand that going from a 32bit to a completely new 64bit application required a complete re-write but I was assuming that a comprehensive tutorials on migrating from 10.X Add-ins to Pro Add-ins would be made available by the vendor as was the case in the past. Perhaps I am not looking in the right places; I do find the ESRI sites rather difficult to navigate. Does anyone know of such tutorial?

ESRI confirms that legacy ArcObjects are not compatible with Pro stating that a complete re-write is needed but the step-by-step how-to seems to be missing. Amusingly, there is a section titled "Migration from 10.x" in the linked github wiki but it says nothing about migrating at all, rather it answers a short list of with a resounding "No"


texts:
I am writing a QGIS Python Plugin. I am trying to write data to a file using

with open("filename.txt", "w") as f:
    f.write("file text.")
QMessageBox.information(None, "Test", "Written to file.")

It runs the above and displays the message but it does not write to the file. Why is it not writing to the file? No exceptions are raised.


texts: How to Open a File Geodatabase with ogr?

texts:
Possible Duplicate:
  How to access feature classes in File Geodatabases with python and gdal?

I have installed:

python-2.7.msi

GDAL-1.9.1.win32-py2.7.exe

import ogr
driver = ogr.GetDriverByName("FileGDB")
print driver

I receive None

Do I need to install something extra to open GDB?


texts: Difference between featureid and featurenumber

texts:
I have created an Atlas and I want to show the page number of each Atlas page. For that, I can create a label and either use

[% @atlas_featureid %] or [% @atlas_featurenumber %]

But what is the difference between the two? Which one is better to use?


texts: Spherical Mercator projection invented by Google and not conformal?

texts:
I had some certitudes on the Google Spherical Mercator also known under EPSG code 900913 and 3857 before looking at this presentation, I was sure that this projection was in the Mercator family (so conformal).

With this input, I'm now not sure about conformal status of Spherical Mercator. I didn't dive enough in mathematical "stuff" here to be able to qualify validity of the doc.

I also want to confirm / infirm this exact projection was never there before Google Maps landing.

Any opinions on those 2 questions?


texts: Need Help Creating Query Layer From Oracle For Past 365 Days

texts:
I have several Query layers that work fine for generic, long term date ranges such as:

This shows all entries that have the previous YEAR in the date field.

AND To_char(ds.calendar_date, 'YYYY') = To_char(( Extract(year FROM add_months(sysdate, -12)) ))

This shows all entries of previous 12 months as well as ALL of previous year

AND To_char(ds.calendar_date, 'YYYY') = To_char(( Extract(year FROM add_months(sysdate,-12)) ))
AND To_char(ds.calendar_date, 'YYYY') = To_char(( Extract(year FROM sysdate) ))

I need something that will show the previous 365 days.

I've tried using AND To_char(ds.calendar_date) >= To_char(sysdate-365) but I get nothing. I've done several iterations with (and without) the To_char operator but nothing does it. The sysdate-365 gives me records from 2004?! I feel like the answer is simple and it will be something along the lines of an ESRI / Oracle / Query layer oddity.


texts: Using ArcGIS spatial statistics tools on average values?

texts:
I'm wondering if anyone can tell me whether using the spatial statistics tools such as Getis Ord Gi* "hot spot" analysis in ArcView is advisable when your input data is an averaged value, e.g. the average income in a census tract. I don't know enough about the mathematics going on under the hood in the spatial statistics toolset to determine whether using averaged census data would introduce error into the analysis.

I'm looking for clusters of high income and low income dissemination areas (slightly smaller than census tracts) within a Census Metropolitan Area. Ideally, I would like to create a series of maps that compare the location of high/low income dissemination areas across 3 census years, but median income values are not provided for one of them, meaning I would have to rely on average income in order to do so. I have run the Getis Ord Gi* tool on average income for these three censuses, as well as median income for the years where it is available and the spatial patterns produced do not seem to differ much. However, this question is still of general interest, as I haven't been able to track down any studies that use Hot Spot Analysis directly on average income for comparison.


texts: How to connect QGIS with Access database

texts:
I am trying to connect Access database with QGIS like I would in ArcView but I get this problem.

Data source name not found and no default driver specified QODBC3: Unable to connect

Hoping someone can help.


texts: How to create TOC icons from a geojson leaflet style

texts:
I have a leaflet map that shows a couple of geojson layers. The style for each layer is represented in a TOC by means of the plugin Grouped Layer Control. The TOC has the following aspect:

The question is if there are tools to generate the icons of the layers for the TOC more or less automatically from leaflet styles like these ones:

function styleMessalo_T100(feature) {
    return {
        weight: 0.52,
        color: '#000000',
        fillColor: '#3da06d',
        opacity: 0,
        fillOpacity: 0.3
    };
};

function styleEstradas(feature) {
    var highway = feature.properties.tags.highway;
    if ( highway && ((highway === 'primary') || (highway === secondary'))) {
        return {
            weight: 0.6,
            color: '#9c0000',
            opacity: 1
        }
    } else {
        return {
            weight: 0.6,
            color: '#6d7672',
            opacity: 1
        }
    }
};

var iconBarragem = L.icon({
    iconUrl: '../markers/barragem.jpg',
    iconSize: [20.0, 20.0],
});

function doPointToLayerBarragem4(feature, latlng) {
    return L.marker(latlng, return {
        icon: iconBarragem,
    });
}

Now, I'm creating similar styles manually in qgis, and using qgis2web plugin to generate a webmap. Then I take the png files from the legend folder, rename it, and paste it into my own project.

To clarify. I'm aware that I can create a more or less pure html+css legend pretty similar to the desired styles, coding it by hand as exposed in this or this or this answers.

But I'm looking for a gui tool or script (preferable for automatizing the process) that given a leaflet style or something like mapbox-simplestyle generates an image file or html+css code that represents the legend.

I also see a suggestion to use geoserver but this seems overkilling.


texts: Calculate raster landcover area in multiple overlapping buffers

texts:
I have 900 1km buffers (polygons) in one layer and a raster layer for forest area (reclassified into forest 1, non forest 0).
I have tried:

Tabulate area - but realised this does not work for overlapping polygons
Convert Raster to Polygon - Intersect and then calculate geometry for area but I'm getting areas equal to the pixel size of the raster (100mx100m= 10,000) and I'm not sure how to convert this back to m2...


texts:
How can I add custom paper sizes to the presets under the Composition tab?

If I save and reuse a template, only the Item properties are saved, not the Composition items, so I have to change the width and height of the paper every time.


texts: What format is this GPS tracking string?

texts:
I am receiving tracking data from a device and I am not sure what is the format/standard of this string in order to properly process the variables, this is a sample string received:

#355488020131816##1#0000#AUT#01#2260012b9cbb47#2112.703000,E,4545.754120,N,0.07,0.00#031213#184625.000##

I am not sure about the variables 2,3,4,5, (##1#0000#AUT#01#) the others seem obvious, position, speed, bearing, date,time.

The device appears to be GT02A


texts:
I'm using QGIS 2.18.16 and following the training manual, and I'm stuck on 6.4.5.

I've followed the example to produce an action that will search wiki, but I receive a syntax error as listed below:

An error occurred during execution of the following code:
from PyQt4.QtCore import QUrl; from PyQt4.QtWebKit import QWebView; myWV = QWebView(None); myWV.load(QUrl(’http://wikipedia.org/wiki/Swellendam’)); myWV.show()

  File "", line 1
    from PyQt4.QtCore import QUrl; from PyQt4.QtWebKit import QWebView; myWV = QWebView(None); myWV.load(QUrl(’http://wikipedia.org/wiki/Swellendam’)); myWV.show()
                                                                                ^
SyntaxError: invalid syntax

The code:  from PyQt4.QtCore import QUrl; from PyQt4.QtWebKit import QWebView; myWV = QWebView(None); myWV.load(QUrl(’http://wikipedia.org/wiki/[% "name" %]’)); myWV.show()

Any thoughts on where I've gone wrong?


texts: GeoWebCache: Is there a way to empty the cache only for a certain zoom level?

texts:
This is a follow-up to this question. In short, I was overambitious in seeding and now my /usr filesystem is full. Take a look at the cache occupancy broken down by zoom level. (Data is OSM South America slice by Geofabrik, imported by osm2pgsql.)

[root@opengeosa opengeo_southamerica]# du -sk EPSG_900913*
4       EPSG_900913_00
8       EPSG_900913_01
16      EPSG_900913_02
36      EPSG_900913_03
236     EPSG_900913_04
832     EPSG_900913_05
3236    EPSG_900913_06
11088   EPSG_900913_07
32828   EPSG_900913_08
110464  EPSG_900913_09
392380  EPSG_900913_10
1487876 EPSG_900913_11
5303356 EPSG_900913_12
20030184        EPSG_900913_13
15776620        EPSG_900913_14
8828    EPSG_900913_15
189778220       EPSG_900913_16
7932    EPSG_900913_17
3764    EPSG_900913_18
492     EPSG_900913_19
200     EPSG_900913_20
352     EPSG_900913_21
[root@opengeosa opengeo_southamerica]#

So, I would like to empy the misbegotten level 16 cache and leave all others. Is this doable? Could I just rm -rf EPSG_900913_16/* or would that be a Very Bad Thing™? Or should I just empty the entire cache and reseed only through level 14?


texts: Fill or inserting new field based on comparison of two other field

texts:
I have a layer Network with two different columns

Route and section

I need a query or function that compare Route rows , and section rows, and create or rather fill the new_field with this condition:

IF(route.row=route.row-1;"newfieldvalue"=1;
IF(section.row=section.row-1;"newfield_value"+1;1))

This Query or function should compare each Coloumn with it selfs that mean:

If the row value of Route is equal to the one above; and if the row value of section couloumn is equal to the one above than insert( @current value of new field + 1) else 1;


texts: Load shapefile in PostgreSQL with intersection rule

texts:
I have one segmentation shapefile vector that I would like to load in a PostgreSQL database. However I need to remove the polygons created in background values of landsat image. Since I have a big amount of features, the performance must be good.

Only the polygon that intersect another polygon can be loaded into the postgres table. This polygon is the same as useful area of image.

Is there any chance to do this operation with spatial restriction when save shapefile in PostgreSQL database?

I need to do it with Python 3.

The first image shows the size of my problem.


texts: Reading Synthetic Aperture Radar raster images using python

texts: Is it possible to read and modify Synthetic Aperture Radar raster images using python? gdal is not able to understand their Type which is Binary (application/octet-stream). One can find sample radar images here.


texts: Switch between test and production database in ArcMap?

texts:
I have two Oracle 11g databases, one for testing (tst.world) and one for production (prd.world). Except when I am trying out new things, they have exactly the same table structure.

In ArcCatalog I have created a connection - let's call it My Connection - to the test database. I use it in a .mdx file to display data from multiple tables in the test base. When the file is opened I get a prompt asking me to enter username and password:

I would like the user to be able to choose between the test and the production database in that prompt, or in some other way choose which database to get the data from. Right now this is not possible since the field "Instance:" is greyed out. So far the only way in which I have been able to accomplish this is to use two separate .mdx files whit exactly the same content except that they use different connections. Clearly that is not a good solution.

Another solution would be to edit the connection settings in ArcCatalog every time I want to switch database. That works for me, but I am afraid it would be considered too complicated and tedious by many of the end users.

Is there a better way to accomplish this task?


texts: Select features that share the same value in a certain column but also have different values in a different column

texts:
I have a dataset of plants, they are many different plant species that are given a numerical code in a column in the attribute table, but they are located in different land types which was also identified in another column , I want to make a selection that highlights the plant (or points) that share the same species ID but are located in different land types. Here is an example:

I want to know if there is a tool that I can use to select the points with OID 2, 3, 7, 8, and 9 since they are grouped by plant species but there is at least one point (i.e. plant) located in different land Types than all the other points that share the same Species ID. I am using ArcMap 10.2.1.
Quick note: there is no map or layer for the Land Type field.


texts: Transferring custom QGIS settings from one computer to another?

texts:
Can I export custom general settings from my QGIS installation to a file so that I can transfer the settings to QGIS on a different computer?

Does QGIS have them already saved to a file located somewhere on hard disk?