texts:
Python Version 3 - Qgis 2.14

1 / what's wrong ?

Lot's of solutions are exposed in differents forum but any work for me.
I need to get path to save documents later and i had this error :

  TypeError: QFileDialog.getOpenFileName(QWidget parent=None, QString
  caption=QString(), QString directory=QString(), QString
  filter=QString(), QString selectedFilter=None, QFileDialog.Options
  options=0): argument 1 has unexpected type 'instance'

Here is my code :

# set event when user click
# param - function with object argument
dlg.resBtn.clicked.connect(self.getFolder(dlg.txt_res))

# function call by click to open dialog
# param - obj : textEditor object update by directory selection
 def getFolder(self,obj):
        ddir= QFileDialog.getExistingDirectory(self, 'Select directory')
        obj.setText(ddir)

2 / what is this error ?

If i give no param to QFileDialog.getExistingDirectory() function, windows to explore is open when code is reading. but i just need to open this windows when user click on push button.

If a close window without select directory, i had this error :

  dlg.resBtn.clicked.connect(self.getFolder(dlg.txt_res))     TypeError:
  connect() slot argument should be a callable or a signal, not
  'NoneType'


texts: Where is this raster offset being introduced in the QGIS processing toolbox?

texts:
I'm using the (awesome!) processing toolbox in QGIS 2.0 to reclassify some NoData values in a raster from -9999 to 0, but at some point in the process there's an offset of about 0.2 degrees being introduced to the output layer.

I'm using the "Reclassify (Simple)" tool (I believe it's part of SAGA), with some pretty basic parameters:

. . . and where the input raster was situated correctly in comparison to this hex grid:

. . . the output raster is bumped southeast by two hex cells, or about 0.2 degrees:

The input origin is -180, 83.623 . . . but the output origin is -179.05,83.05. Normally I would assume this is the result of a faulty WGS84 conversion, but the input, output and project are all EPSG:4326.

Where is this origin offset coming from, and is it worth filing a SAGA ticket or is it just something I overlooked?

Update
The process log indicates a transformation is performed on the input origin:

Transformation x' = -179.999989 + x * 0.100000 + y * 0.000000
Transformation y' = 83.623125 + x * 0.000000 + y * -0.100000

and then a few lines later the new parameters show up, now incorrect:

Grid system: 0.1; 3600x 1735y; -179.949989x -89.826875y


texts: I have a point shape (order number, value) which I want to convert in poliline by order of points from 1,2,3,.....x, i try in QGIS plugin Points2One but I did not succeed, can someone help me with any sugestion?


texts: Populating the drop down of GPValueTable

texts:
In a Python Toolbox the user has to select a city from a list of 50 world cities. More than one city can be selected.

The best UI control for that would be Autocomplete / Autosuggest, but since this does not exist I thought of a GPValueTable with a single column of type String. This works well when the user types in the name of the city and then clicks the plus button:

And it's possible to validate it against the list of cities using updateParameters. However, I would like the drop down list to be populated too:

The question is how can the drop down list be populated with a custom list of strings?

The code in the getParameterInfo function:

param0 = arcpy.Parameter(
displayName="Cities",
name="in_features",
datatype="GPValueTable",
parameterType="Required",
direction="Input")
param0.columns = [['String', 'Name']]

Update:

There is the option of using the out-of-the-box multi select control multiValue=True creates many checkboxes and the user has to scroll. Typing in the keyboard only finds items by the first letter.


texts: Using WCS layer in QGIS for processing gives FAILURE: Missing destination

texts:
I have a GeoServer serving raster dem data out via WCS and can add the layer easily to QGIS (3.2) but when I try to use it as the input for an operation (e.g. HillShade) I get an error.

Processing algorithm…
Algorithm 'Hillshade' starting…
Input parameters:
{ 'ALTITUDE' : 45, 'AZIMUTH' : 315, 'BAND' : 1, 'COMBINED' : False, 'COMPUTE_EDGES' : False, 'INPUT' : 'dpiMode=7&identifier=training:panorama&url=http://localhost:8080/geoserver/ows?version%3D1.1.1', 'MULTIDIRECTIONAL' : False, 'OPTIONS' : '', 'OUTPUT' : 'C:/Users/ishare/AppData/Local/Temp/2/processing_868cc25180f24d9f8e193bc0cfedb6d3/7cf48f64d8b54ee694c3ee93ed17e903/OUTPUT.tif', 'SCALE' : 1, 'ZEVENBERGEN' : False, 'Z_FACTOR' : 5 }

GDAL command:
gdaldem hillshade dpiMode=7&identifier=training:panorama&url=http://localhost:8080/geoserver/ows?version%3D1.1.1 C:/Users/ishare/AppData/Local/Temp/2/processing_868cc25180f24d9f8e193bc0cfedb6d3/7cf48f64d8b54ee694c3ee93ed17e903/OUTPUT.tif -of GTiff -b 1 -z 5.0 -s 1.0 -az 315.0 -alt 45.0
GDAL command output:

FAILURE: Missing destination.

Clearly the missing destination is a red herring, and it is failing to understand the WCS layer as an input. But GDAL is built with WCS support included so shouldn't this work?

Do I need to do something else or start working out where to hack the QGIS code?


texts: Missing optional dependency 'tables' in QGIS

texts:
When we run a Python script as an action in QGIS, we get an error as shown below.
Traceback (most recent call last):
  File "", line 17, in
  File "C:\PROGRA~1\QGIS3~1.16\apps\Python37\lib\site-packages\pandas\io\pytables.py", line 397, in read_hdf
    store = HDFStore(path_or_buf, mode=mode, errors=errors, **kwargs)
  File "C:\PROGRA~1\QGIS3~1.16\apps\Python37\lib\site-packages\pandas\io\pytables.py", line 518, in __init__
    tables = import_optional_dependency("tables")
  File "C:\PROGRA~1\QGIS3~1.16\apps\Python37\lib\site-packages\pandas\compat\_optional.py", line 92, in import_optional_dependency
    raise ImportError(msg) from None
ImportError: Missing optional dependency 'tables'.  Use pip or conda to install tables.

We installed the tables and we have tried to add the paths that QGIS is using from Environment Variables like in this question (https://stackoverflow.com/questions/58479748/missing-optional-dependency-tables-in-pandas-to-hdf) but nothing has changed. Unfortunately, I do not know how to use conda while creating an action in QGIS.
We have already installed the required packages in OSGeo4W Shell and did not get an error/warning message.
C:\>pip install tables
Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: tables in c:\users\semiha - pc\appdata\roaming\py
thon\python37\site-packages (3.6.1)
Requirement already satisfied: numexpr>=2.6.2 in c:\users\semiha - pc\appdata\ro
aming\python\python37\site-packages (from tables) (2.7.3)
Requirement already satisfied: numpy>=1.9.3 in c:\progra~1\qgis3~1.16\apps\pytho
n37\lib\site-packages (from tables) (1.19.1)

What can we do to solve this problem?


texts: Setting data source for certain layers in .mxd using ArcPy

texts:
I have an .mxd that I use as a template in ArcMap 10.8. I currently manually update the data source for certain layers in the MXD. I am trying to create an arcpy script to do this for me but am stuck at the part of how to change the data source for these layer(s).
I believe that I need to first determine the source of the layers I am changing and then use that to change the other layers with that same source. I need to do this as sometimes my template sources get changed by accident so I cannot put a set path in my script. I am just unsure of how to determine the source of a layer and how to change it and the source of the other layers in my MXD that have the same found source.
Here is a picture of the layers I change.

The application Tenure layer is a polygon fc and tenure road application layer, Pofc and PofT layers (all have the same source) are line fc's. If someone could just help me with the code to find the source of one of the layers, I can figure out the loop to select which layer to find the source from. And I could also use some help on how to replace the source of other layers with that source. I believe I have to use:
Layer.replaceDataSource(workspace_path, workspace_type, dataset_name, {validate})

If yes, I am unsure of how to use it.


texts: How to isolate certain raster values using raster calculator?

texts:
I have a one band grayscale image. It has values 0-63 within the band. I need to isolate some of those values. Is there a way to do this in the raster calculator?

QGIS 2.2//Windows 7


texts: Using R to count number of pixels in raster file for which there are at least one point object present

texts:
I am trying to count the number of pixels in a raster file for which there are at least one point object present in R. I think this would be easy but I haven't been able to figure it out.
I am new to spatial analyses and finding the proper terminology to describe what I am hoping to achieve is challenging.
A Visual Representation of the Problem
Here we have a hypothetical raster with a multi-point vector overlaid. There are 10 pixels in the raster, 3 of which have points associated with them.

Ideal Outcome
I want to write some code in R that, given this hypothetical scenario, would return a value of 3 indicating that 3 of these pixels in the raster file have a point object present.
Do you have any tips or tricks?


texts: How can I display a symbol to a specific attibute in QGIS

texts:
I have a point shapefile, where I categorized my waypoints with attributes (schools, toilets, etc) and I would like to display certain attributes with specific symbols.

How can this be done?


texts: Visualizing Strahler-Order bigger than a specific value

texts:
I try to visualize the Strahler-Order bigger than 8 using raster calculator.

I changed the style of the raster by changing the render type to singleband pseudocolor and, in the colors selector, I did not find boolean sentences. So, I created a new color map and classified it into two classes between 0 and 1 values.

However, it shows me just one small part of my raster (see Figure). The size of the data is as big as the original file. So, I assume that the calculation was properly carried out.

Any suggestions??


texts: Identifying which tracks are too close together?

texts:
I have a point shape (yield logging data) and I want to find the tracks that are too close to another track. Usually the combine drives always with the same distance (e.g. 10m). I want to find the tracks that are closer then 10m (e.g. 5m). Could this work with the voronoi polygons?

The picture below showes the problem. The red market distance between the points is ok. The blue marked point distance is not ok. These points I want to find/select.


texts: How to store large amounts of text in an SDE table

texts:
I need to store a lot of text in a field in an SDE registered business table (an ITable, not a featureclass) - let's say 10,000 words roughly.  I am using Oracle 10g and ArcSDE 9.1.

I am programming with ArcObjects 9.2 sp6.

It looks like the only option might be to use a BLOB as a CLOB is not usable with 9.1 (shame).
Can you use a BLOB in this way and how would you do it in ArcObjects code?


texts:
Is the Density View for CartoDB no longer working? I have looked at three datasets on multiple accounts, and none of them are working correctly. They show only 1 bucket color, no matter the density. The CSS is only generating one grouping for the point density, so obviously it is not generating the correct CSS. I can manually manipulate it (the CSS), but requires a lot of fine tuning to find the right grouping!

If anyone knows something I could be doing wrong, please tell me!


texts:
I have a SRTM-3 file and I want to use it to make a DEM of the related world region. I'm using MapInfo Professionale V8 and I'm fairly new to that software,

so if there is any tutorial that shows how to use the SRTM file to get a DEM in order to make an elevation map it would be great.


texts:
Will the ArcGIS frequency analysis tool work on a shapefile's attribute table or does it only work on a discrete/standalone  table loaded in a mapping session?

Have loaded a geolocated address book into map and am trying to get a count by zip code (one to many), so I might join count to zip code shapefile and symbolize by "quantity".


texts: Can I use EPSG:25832 and EPSG:32632 interchangably

texts:
up until now we have been using EPSG:25832 (UTM Zone 32N) as our standard Coordinate system. Recently however after we aquired a new project area and consequently new software, I have to display data with EPSG:32632. Looking at the definitions, it seems the only difference is between the extents and the Datum. EPSG:25832 is based on the "European_Terrestrial_Reference_System_1989" using the spheroid "GRS 1980" and EPSG:32632 uses the datum "WGS_1984" with the spheroid of the same name - This means the Semi-major axis are the same but the Semi-minor axis has a difference of 2 meters. All other parameters are the same.

So now the question....can I use them interchangeably on a regional scale? meaning I will not be using them to control satellite systems or aircraft flight paths.


texts: Publishing Map as a Feature Service Ends up as a Map Server

texts:
Following ESRI documentation at Here

I am trying to publish a Map as Feature Server by enabling the Feature Access on Capability bur when I publish the Map the result looks like

As you can see the URL ends with /MapServer while I need to get a /FeatureServer like

Can you please let me know what I am doing wrong?


texts:
Can you suggest how a JasperReport can be generated by an user accessing an point (say a wireless tower in map) on a active layer done in OpenLayers 3. Here the architecture uses OpenLayers 3 and is linked with GeoServer and PostGIS.

The jasper report that needs to be launched for the selected wireless tower or point, needs to pick some spatial coordinates (location, intersections etc) and non-spatial information/data (installation date of tower, some runtime statistics info of tower etc) from PostGIS database, for the respective telecom tower.

Do share your views.


texts: DEM and vector data won't match using QGIS

texts:
I have just started using QGIS for my research at university. I have a DEM comprising Italian regions of Lazio, Tuscany and Umbria. I then have a set of vector data from surveys carried out in these territories. When I upload the DEM and survey data, I expect them to be on top of one another, but they are not. I tried everything but they just won't match! If I put a base map from Google Earth or Maps, the survey data matches the Italian regions they refer to, while the dem data ends up corresponding to Spain for some reason. I changed the coordinate references, to Monte Mario Italy 1, but nothing really changes: actually, the survey data just disappears from the map! It's probably something really stupid that I am messing up, but I tried every possible combination and I am getting quite frustrated!

Does anyone know what I can do to fix this?


texts:
I have two Landsat images for the same study area. One of them in 1975 (Landsat 2 MSS with 60m resolution) and the other in 2014 (Landsat8 OLI with 30m resolution). I want to do change detection between them. As you know in order to do change analysis on images with different resolutions, I should make them all the same resolution so they match. Since Landsat 8 is higher resolution from Landsat 2. So, I need to do resample for scene of Landsat 8 to to match the MSS data and to compute the change detection.

So please, currently I am working on ENVI 5.1 software. How can I resample the Landsat 8 and do change detection between two satellite images in this software?


texts: Clipping Shapefile of Circle features using Raster in ArcGIS Desktop?

texts:
I want to keep all the features in the shapefile that lie in the raster region and reject others. How do i do this?

Data to be Clipped: Is A Shapefile of 1 mile circles.

Data I want it clipped with : is a Raster (GeoTiff or ADF).

Now I want to find out all the Circles from shapefile that overlap on the Raster.

I do not want to modify the raster here. I want those features (circles) from the Shapefile removed which do not overlap with the raster.

I have tried the tool mentioned by Dano, but tools expects a feature input in the second field, and not a raster directly.


texts: how to aggregate raster values for each grid in a fishnet-like map

texts:
here is my problem:

I have a raster layer containing the elevation info. and another gridded layer (shapefile type). What I want to do is to calculate the mean elevation within each grid to have a single representative value.

Does anyone know how to handle this? Thx in advance for any suggestion .


texts: Polygon interpolation neighbour average?

texts:
I am trying to estimate empty values based on a polygon layer we received from another study. Is there a way to calculate an average value for polygon A based on neighbours B,C, etc that border it?

I thought about converting polygons to points and interpolating, but thought there may be another method of doing this.

Photo attached for reference


texts:
I am working on project In which I am supposed to make a web app which will be asking users to upload their inputs as "Rinex file" of Rover and Base Station and It will give more precise data after differential correction.

So I just want to know the algorithm of post processing kinematics.


texts: ArcGIS API for Silverlight & Microsoft Active Directory Security

texts:
I am trying to get a Silverlight app working with an Esri map service, secured with Active Directory authentication.

I found an example online of how to get the app working when using User Store authentication at https://developers.arcgis.com/silverlight/sample-code/start.htm#IdentityManagerServices.

Does anyone know how I can get the app working correctly when the user has a map service with Active Directory authentication?


texts: Using North Pole Orthogonal projection in QGIS?

texts:
I was looking for the North Pole Orthogonal projection in QGIS but I haven't found it yet.

Is there any way to add this projection in CRS?


texts: Setting defaults after a value is chosen

texts:
I am developing a way for field workers to use arc collector on tablets.

I ran into a problem during development. Here is the way its set up: three workers create points in arc collector on three separate tablets, and in the attribute table for these points is a field to input their initials to show who created the point. I have created a drop-down of their initials to make it easier, but they are saying it would save time to have the default set to their initials instead of having to use the drop down for every point. I do not know how to set a default for three different initials for three different tablets. Is it possible to set a default in a domain after a value is chosen instead of before, even if its just in ArcMap? Is there another way of getting around this?


texts: ArcGIS Online related table not displaying correctly in popup

texts: I have an ArcGIS Online Map composed of three spatial feature classes and a table.  All three feature layers are related to the table and all the layers and table are housed within the same SDE.  I have enabled the popups for each layer using fields from the related table in the popup.  However, instead of showing the correct values from the related table I only see 1s and 0s as values.  If the value in the related table is NULL already the popup shows a 0 and if it is not NULL it shows a 1.  How can I resolve this issue so I can see my actual values in the popup?[


texts: Looking for average temperature data for Afghanistan

texts: Hey all I am looking for a average temperature data for Afghanistan monthly over a year, precipitation would work too. In a shapefile, if anyone has any or knows anywhere that I could access one please let me know.


texts: Join by attribute / spatialite with SQL / left outer join in QGIS

texts:
I have a layer with polygons (let's call it "woods") that is intersected by "n" line-objects ("roads"). Every road has an unique attribute roadid.

I need to add ALL the roadid's (not only the first QGis finds) to the woods-polygons-layer they intersect for latter use. Preferably all the roadid's should be in one new attribute column of the "woods", divided by e.g. ",".

If there was only one road in every wood i could use the "join attribute by position" tool to get roadid's on the polygons. The wanted attribute is an unique string so it doesn't help to sum/average/min/max the fields, and there is no option to bind them together as strings divided by ",".

Edit:
There is no solution without the help of plugins or external programms by now (see comment by Matthias Kuhn below).

(QGIS 2.10.1-Pisa)


texts:
I wrote this pyhton script to clip a raster with a shapefile:

import ogr
import subprocess

inraster = '/Users/sigc2sige/Desktop/gdal/rast.tif'
inshape = '/Users/sigc2sige/Desktop/gdal/poly.shp'

ds = ogr.Open(inshape)
lyr = ds.GetLayer(0)

lyr.ResetReading()
ft = lyr.GetNextFeature()

while ft:

       com = ft.GetFieldAsString('com')

       outraster = inraster.replace('.tif', '_%s.tif' % com.replace(' ', '_'))
       subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape,
                 '-crop_to_cutline', '-cwhere', "'com'='%s'" % com])

      ft = lyr.GetNextFeature()

ds = None

When I execute the code I get this error : FileNotFoundError: [Errno 2] No such file or directory: 'gdalwarp': 'gdalwarp'

I have installed gdal but I don't know why the gdalwarp is not recognized. I'm running on Mac.

Best


texts: List of lat-long, need values in a raster

texts:
Apologies in advance as I am not a GIS specialist by any means. I have a set of 1 million points and I'm trying to find their values in a geotiff raster. I've tried various versions of this answer involving affine transformations:
https://gis.stackexchange.com/a/221471/143163

def retrieve_pixel_value(geo_coord, data_source):
    """Return floating-point value that corresponds to given point."""
    x, y = geo_coord[0], geo_coord[1]
    forward_transform =  \
        affine.Affine.from_gdal(*data_source.GetGeoTransform())
    reverse_transform = ~forward_transform
    px, py = reverse_transform * (x, y)
    px, py = int(px + 0.5), int(py + 0.5)
    pixel_coord = px, py

    data_array = np.array(data_source.GetRasterBand(1).ReadAsArray())
    return data_array[pixel_coord[0]][pixel_coord[1]]

This isn't ideal since it's point by point querying but it's better than nothing. The problem I'm having, however, is that it's not returning the correct values.

As a sanity check, I used some WRF data that had lat/long layers and queried various points and the resultant lat/long that is supposed to be the nearest coordinates in the WRF data are very far off from where they should be. e.g. inputting 33.77864,-117.33142 returns 38.72556, -115.75209 (the range of this layer is 32.597065 to 39.3944,-121.413025 to -113.04607, so there should be much closer matches). Furthermore, switching lat long as inputs doesn't drastically change the return value (when switched, returns 34.820377,-120.55661). Like I said, not an expert, but to me this seems like I'm using the wrong coordinate system as inputs. Does anyone know a way to convert lat long to the appropriate coordinates to find values in a raster?

I realize this isn't the most efficient way to do a list query on a big db, given the original problem of finding raster values for 1 million points, is there a more GIS-ish way to do this?


texts: Displaying mosaic dataset footprints with ArcGIS Image Service

texts:
Is there a way to display the footprints created through a mosaic dataset with an Image Service?

I have a mosaic dataset that created footprints with associated data that can be viewed in ArcCatalog. However, when shared as an image service, only the image is displayed.  You can open properties and make the footprints viewable, but this doesn't show up in the table of contents for selection or detailed functions/processes.

Does anyone know how to make the footprints a selectable layer in an Image Service, or is this the wrong way to go about doing this?


texts: Creating Regular points every meter at county level using ArcGIS Desktop?

texts: I have a need to create regular points with a coordinate every meter in a county. I've tried several methods including Hawthes Tools in ArcMap9.3, ET GeoWizards, and the standard Toolbox. I have encountered reoccurring instances where all of the processes have froze up.I am assuming that these programs can't handle that data generation. I tried this in QGIS as well. I need help ASAP.


texts: How to fix JSON error in Silverlight Application?

texts: I am trying to run a simple Silverlight application based on code provided by ESRI at this link (Silverlight Samples).  My app works fine if I just use the TiledLayer but if I add the Dynamic layer and set it to a Wells MapService I have running on my server, I get an error stating that the WellsLayer has a SerializationException and is not a valid JSON primitive and this error can occur when extraneous data is present after JSON data.  What does that mean?  I did see in the forums that someone suggested that a clientaccesspolicy may not be in place at the wwwroot folder.  I have downloaded and placed that file there but no luck.. Thoughts?


texts: Simplified lines disappear in smaller zoom levels in ArcGIS JavaScript v.4.14 map

texts:
I have a GraphicsLayer on ArcGIS js map with many lines of varying lengths (1 to 200m).
However, in smaller zoom levels the short lines disappear. For example, in the images below the red lines are visible in the larger zoom level and then they are no longer visible when zoomed out:

ArcGIS zoomed in

ArcGIS zoomed out

Having the same layer on Leaflet does not have this effect and the short lines appear like points:

Leaflet zoomed out

Here is the code:

    this.simpleLineSymbolLeft = new SimpleLineSymbol({
        type: "simple-line",
        color: [247, 45, 45],
        width: 5,
        style: "solid"
    });

let graphic = new Graphic({
            geometry: polyline,
            symbol: this.simpleLineSymbolLeft,
            attributes: properties
        })

I assume this is a simplification effect, based on which very short line, instead of being collapsed into points, they are turned off. Is there a way to tweak this so that the short lines will not disappear?


texts: Is the ServerAction Method Required in a ImapServerCommandAction?

texts:
I am using an ESRI toolbar control buddied with a map control.

I am using a command.  Below is a simple example.

My question is to clarify - having the ServerAction method is a requirement when you are using the toolbar control and server side action?

public class CustomExtent : IMapServerCommandAction //implements the IServerAction interface
{
    public CustomExtent()
    {
    }

#region IServerAction Members

public void ServerAction(ToolbarItemInfo info)
{
   Map map = info.BuddyControls[0] as Map;
   map.Extent = new Envelope(-90, 10, -70, 5);
}


texts:
I have a vector layer in QGIS with 30000 points. I want to delete those with values lower than a specific value. I could not use filter in an edited layer, and I could not filter the layer after editing it.

What can I do?


texts: Add a Mapinfo TabFile in Geotools goes wrong [Automatically Converts my Integer values to Byte]

texts:
we have read Mapinfo Tab file in Geotools using JniOGRDataStoreFactory, but when we execute filter statement on feature collection it give some vague values so we decide to investigate it and found that it create some attributes with byte data type that originally are integer.

Console Image from Eclipse After Reading using JniOGRDataStoreFactory

Mapinfo tools screen shot.

now my problem is that how can we read that attribute values correctly, Please advice.


texts:
I have setup a draggable marker on my map with Leaflet API. Its working fine. But the problem is when I am trying to drag it to the border area the map is not moving. I want a panning map with Marker dragging. Please let me also know how can able to control the map panning when dragging the marker.

Thanks


texts: arcpy.Project_management is projecting points but not polygons?

texts:
In my current project, I am grabbing data from a web service and want to change the projection before inserting into a production db. For some reason when I am projecting the data, some of it is being lost. The point feature class is being projected just fine, but the polygon feature class created after calling arcpy.Project_management has no data in it. I've created sample code of the issue and hope somebody can shed some light.

import arcpy
import os

# scratch gdb and 2 feature classes (these are WGS84)
scratch_gdb = os.path.join("C:\\", "temporary", "scratch.gdb")
original_points = os.path.join(scratch_gdb, "original_points")
original_polys = os.path.join(scratch_gdb, "original_polys")

# remove rows from our feature classes
arcpy.DeleteRows_management(original_points)
arcpy.DeleteRows_management(original_polys)

# add a point to the point feature class
point = arcpy.Point(-80.247928, 33.464786)
cursor = arcpy.da.InsertCursor(original_points, ['SHAPE@XY'])
cursor.insertRow([point])

# add a polygon to the polygon feature class
polyarray = arcpy.Array([arcpy.Point(-80.254769, 33.467541),
                         arcpy.Point(-80.239318, 33.468623),
                         arcpy.Point(-80.236017, 33.456030),
                         arcpy.Point(-80.264083, 33.451404)])
polygon = arcpy.Polygon(polyarray)
cursor = arcpy.da.InsertCursor(original_polys, ['SHAPE@'])
cursor.insertRow([polygon])

# at this point, both original_points and original_polys have data

# create path for 2 temporary feature classes
projected_points = os.path.join(scratch_gdb, "projected_points")
projected_polys = os.path.join(scratch_gdb, "projected_polys")

# delete these feature classes if they exist (from previously running the script)
if(arcpy.Exists(projected_points)):
    arcpy.Delete_management(projected_points)
if(arcpy.Exists(projected_polys)):
    arcpy.Delete_management(projected_polys)

# project our 2 old feature classes into our 2 new feature classes
sr = arcpy.SpatialReference(2273)
arcpy.Project_management(original_points, projected_points, sr)
arcpy.Project_management(original_polys, projected_polys, sr)

# projected_points has a point, but projected_polys has no rows ?


texts: ArcMap, equalize the value of xll corners and yll corners in two different rasters

texts:
I have DEM and landcover rasters, I have different values between my DEM and land cover map for xllcorner and yllcorner.
And how can I fix this?


texts:
I recently received a PostGIS instance and am trying to figure out how much of my current processing workflow (in proprietary tools) can be ported over to PostGIS.

One of the most frequent tasks I encounter is validating/repairing geometry.

Here's what I've done and here's where I'm stuck:

I load some data (multipolygon) to postGIS.  Then I use the following statement to identify and copy invalid geometries to a different table:

create table broken_polygons as
select * from polygons
where not st_isvalid(geom);

Then, I get rid of the offending records:

delete from polygons
where not st_isvalid(geom);

This gives me a chance to take a closer look at the bad geometry; sometimes there are only one or two, sometimes there are hundreds.  But when I try to do this:

select *, st_makevalid(geom) as fixed_polygons
from "broken_polygons"

It gives me a result, but when I load up fixed_polygons in QGIS and check the geometry validity (GEOS), all of them error out: nothing has been modified.

In this particular case, all of the errors are self-intersections.

I am using PostGIS 2.3.2 on postgres 9.6.2

I am not a DBA, just someone trying to get things done with the tools in front of me; so I seek examples of precise syntax to get ST_MakeValid to work.


texts: Get Data for ArcPad extension and python

texts: Is it possible to call the 'Get Data for ArcPad' extension using python?  If so would someone be willing to share their code?


texts: How to compare outputs of erosion from ArcSWAT and RUSLE from ArcGIS?

texts:
I calculated water erosion in ArcSWAT and I also calculated erosion by RUSLE

http://35.8.121.139/rusle/

equation in ArcGIS. Now I need to compare these two (from ArcSWAT I used SYLD- assumed sediment routing in t/ha/y and in rusle erosion also in t/ha/y), I used zonal statistics for rusle, but my results are wrong and I think I am using the wrong instructions on how to do it, does anybody know how can I compare these two?


texts: Reclassify raster based on cell location and value

texts:
I am working in ArcGIS 10.5 (with the spatial analyst extension); I'm relatively inexperienced but I understand what I am trying to do, just need some help on how to accomplish it.

I am working with MODIS (raster) data. I am trying to reclassify (group) cells in a raster to have the same value if 1) in the same 3x3 neighborhood AND 2) if the difference of their values is less than 6. See picture below for example of what I am trying to do.


texts: Creating an INFO table using ITable in Arc (Java)

texts:
I'm trying to set up an info table in an addin for ArcGIS 10.  I'm using Java and Eclipse IDE.  The code I have is here:

public class Tables {

String database = "C:\\GIS";
String nameOfTable = "getOn";
public ITable createInfoTable(String database, String nameOfTable) throws IOException  {

        IPropertySet propertySet = new PropertySet();
        propertySet.setProperty("DATABASE", database);

        IWorkspaceFactory workspaceFactory = new ArcInfoWorkspaceFactory();

        IArcInfoWorkspace workspace = (IArcInfoWorkspace)workspaceFactory.open( propertySet, 0);

        //Fails here (or above actually)  Check out Factory.open

        IArcInfoItems arcInfoItems = new ArcInfoItems();
        ITable  myTable = workspace.createInfoTable("myTable", arcInfoItems);

        return myTable;

I've adapted this code from C# as there are no examples on the Arc pages for Java.  You can see with the comment where it is failing.  If you can see what's going wrong please let me know.  And I will repost if I crack it.  (see here for original C# example: http://resources.esri.com/help/9.3/ArcGISEngine/ArcObjects/esriDataSourcesFile/IArcInfoWorkspace.CreateInfoTable_Example.htm)


texts: I want to run I-Tree "Canopy" and am looking for an official Esri shapefile that defines the city of Warren Michigan's jurisdictional boundaries.


texts: How to return more features from GetFeatureInfo?

texts:
I configured WFS service with "Maximum number of features for preview" to 50 but when i try to get features my request returns only one value, for example

{"type":"FeatureCollection","totalFeatures":"unknown","features":[{"type":"Feature","id":"","geometry":null,"properties":{"GRAY_INDEX":-9999}}],"crs":null}

and with &feature_count parameter geoserver accepts only 50 as value but also adding this the ouput is the same the above. Instead I want an output as:

{"type":"FeatureCollection","totalFeatures":"unknown","features":[{"type":"Feature","id":"","geometry":null,"properties":{"GRAY_INDEX":-9999}}, {"type":"Feature","id":"","geometry":null,"properties":{"GRAY_INDEX":4}}],"crs":null}

EDIT1:

In my Geoserver I have multiple time spaced rasters for each layer. When I try the GetFeatureInfo request above I specify also a time range as &time=[from]/[to]. Because the bounding box of these rasters are overlapped, I want that the WMS request returns a list of pixel features when I click on the overlapped area.

EDIT2:

The request is:

geoserver/ows?service=wms&request=GetFeatureInfo&info_format=application/json&feature_count=50&exceptions=application/json&bbox=-4.226473,39.039063,16.516447,79.980469&layers=id210s013016r000-europe&query_layers=id210s013016r000-europe&width=381&height=752&crs=EPSG:4258&time=2016-03-03T23:00:00.000Z/2016-03-04T20:00:00.000Z&i=338&j=134


texts:
I'm trying to implement a Horizontal Slider using 'dijit/form/HorizontalSlider'. I successfully got the Slider working but the slider comes with a red color.

How can i change the color to white or any color that i want?

Below is the code snippet that im using.

 createSlider= function() {
      sliderOpacity = new HorizontalSlider({
          name: "slider",
          value: 60,
          minimum: 0,
          maximum: 100,
          intermediateChanges: true
      }, "sliderOpacity");
      sliderOpacity.startup();
      dojo.connect(dijit.byId('sliderOpacity'), 'onChange', changeOpacity);
    },


texts: How to compute areas where the value in one layer is smaller than in another layer?

texts:
I have two vector layers, one contains actual values (in this case % cover of native species), the other contains target vales of the same attribute. the polygons are not always of the same geometry.

I want to produce a resulting layer of polygons where the actual values are below the target values.

i.e. some sort of polygon spatial overlap where 'layer1', 'actual cover' < 'layer2', 'target cover'.


texts: Black box around raster layer retreieved from Mapbox in ArcGIS online

texts:
I manage a mapping system for the Pine Ridge Reservation in South Dakota and am trying to retrieve a raster layer but I keep getting this expanding black background/box.

When I open the mbtiles file in Qgis it looks fine but when I retrieve the file through ArcGIS either as a WMTS layer OR the Mapbox instructions for pulling layers into ArcGIS online of  "Add > Add Layer from Web > A Tile Layer" there's black box around it that increases in size as you zoom-out.

See image.

Here's some relevant details about my data:

The original file was converted from Geotiff to Mbtiles format in TileOven (active Linux fork of TileMill) and then uploaded to Mapbox and styled in map box studio.
This is not a no-data border in the original geotiff. The original geotiff has only a tiny border that I'm not concerned about concealing.
The original geotiff was created using GDAL processing tools in Qgis
To convert the file from geotiff to Mbtiles in TileMill I used these instructions https://www.youtube.com/watch?v=05Et84xIgls
I made certain I set extends in TileMill and ArcGIS online.
I also tried uploading the Geotiff directly to ArcGIS online and the same thing happens.

I'm using TileOven, QGIS 2.18 on Kubuntu 18.04 with KDE Plasma 5.12.6 Software

Here's the ArcGIS Online integration URL provided by MapBox https://api.mapbox.com/styles/v1/dbartecchi/cjlgwhz4704j02rldjapn5wzc/tiles/256/{level}/{col}/{row}@2x?access_token=pk.eyJ1IjoiZGJhcnRlY2NoaSIsImEiOiJPV0cyMW1NIn0.D5p-ioo0SzDujMVD77ZQIg

Here's the Mbtiles file opened in QGIS 2.18 looking exactly how it should look.


texts:
I have downloaded Land Cover 1990 data from Digimap and clipped to a watershed I have already delineated.

My problem is that my land covers are just numbers and I cannot seem to get labels for the land covers.


texts: OpenLayers - LayerSwitcher MattWalker not working with fetch()

texts:
I have a problem with Matt Walker's LayerSwitcher on my map construct. I use fetch () to get a WFS. Unfortunately, the LayerSwitcher does not respond to the vector and nothing is displayed in the switcher.
Has anyone had experience with this switcher?

fetch('[WFS-URL]').then(function(response) {
    return response.json().then(function(polygon) {...


texts: Sentinal NDVI Time series Chart and CLassification

texts:
I am  working  on NDVI time series analysis using Sentinal data  and also  classifying the data it will give some error and not working and how can i get the chart of three different features in single chart

link script

// Function to mask clouds using the Sentinel-2 QA band
/// * @param {ee.Image} image Sentinel-2 image
/// * @return {ee.Image} cloud masked Sentinel-2 image
function maskS2clouds(image) {
  var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}
var S2 = ee.ImageCollection('COPERNICUS/S2')

//filter start and end date
.filterDate('2018-09-01', '2019-02-10')

//filter according to drawn boundary
.filterBounds(geometry)
 // Pre-filter to get less cloudy granules.
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
.map(maskS2clouds);
// Function to calculate and add an NDVI band
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']));
};
print(S2);
// Add NDVI band to image collection
var S2 = S2.map(addNDVI);
// Extract NDVI band and create NDVI median composite image
var NDVI = S2.select(['nd']);
print(NDVI);
var NDVImed = NDVI.median();

// Create palettes for display of NDVI
//var ndvi_pal = ['#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b',
//'#a6d96a'];
var ndviParams = {min: -1, max: 1, palette: ['blue','white','green' ]};
Map.addLayer(NDVImed, ndviParams, 'NDVI image');

////////////////////////////////////////////////////////////////////////////////////
// Create a time series chart.
var plotNDVI = ui.Chart.image.seriesByRegion(NDVI, geometry,ee.Reducer.mean(),
'nd',500,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'NDVI short-term time series',
                hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'}});

//////////////////////////////////////////////////////////////////////////////////////////////

// Display.
print(plotNDVI);

// Display NDVI results on map
Map.addLayer(NDVI);

var newfc = wheat.merge(must).merge(beng);
// print(newfc);

var training = NDVI.select(NDVI).sampleRegions({
  collection: newfc,
  properties: ['id'],
  scale: 30
});
// print(training)

var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: 'id',
  inputProperties: B1
});
print(classifier.explain());

var classified = image.select(bands).classify(classifier);
Map.addLayer(classified, {min: 0, max: 2, palette: ['red', 'blue', 'green']});


texts: Earth Engine: Error: Image.select: Pattern 'test' did not match any bands

texts:
I am trying to combine two images one containing various bands of the sentinel 2 mission and the other one is created from a feature collection containing two types of polygons that have the values "test" of 0 and 1 as can be seen here:

The combined image should then be exported following the tutorial here (https://colab.research.google.com/github/google/earthengine-api/blob/master/python/examples/ipynb/UNET_regression_demo.ipynb#scrollTo=FyRpvwENxE-A)

Here is my code:

RESPONSE = 'test'

trainingMap = ee.FeatureCollection('users/xxx/xxx')
trainingMapReduced = trainingMap.reduceToImage([RESPONSE], 'max');

# combine the images
featureStack = ee.Image.cat([
  image.select(BANDS),
  trainingMapReduced.select([RESPONSE])
]).float()

# consider a kernel of KERNEL_SIZE x KERNEL_SIZE around each pixel
list = ee.List.repeat(1, KERNEL_SIZE)
lists = ee.List.repeat(list, KERNEL_SIZE)
kernel = ee.Kernel.fixed(KERNEL_SIZE, KERNEL_SIZE, lists)

arrays = featureStack.neighborhoodToArray(kernel)

# load training and evaluation areas

trainingPolys = ee.FeatureCollection('users/xxx/trainingRegion')
evalPolys = ee.FeatureCollection('users/xxx/testingRegion')

trainingPolysList = trainingPolys.toList(trainingPolys.size())
evalPolysList = trainingPolys.toList(trainingPolys.size())

# These numbers determined experimentally.
n = 100 # Number of shards in each polygon.
N = 2000 # Total sample size in each polygon.

# Export all the training data (in many pieces), with one task
# per geometry.
for g in range(trainingPolys.size().getInfo()):
  geomSample = ee.FeatureCollection([])
  for i in range(n):
    sample = arrays.sample(
      region = ee.Feature(trainingPolysList.get(g)).geometry(),
      scale = 30,
      numPixels = N / n, # Size of the shard.
      seed = i,
      tileScale = 8
    )
    geomSample = geomSample.merge(sample)

  desc = TRAINING_BASE + '_g' + str(g)
  task = ee.batch.Export.table.toDrive(
    collection = geomSample,
    description = desc,
    folder = BUCKET,
    fileNamePrefix = FOLDER + '_' + desc,
    fileFormat = 'TFRecord',
    selectors = BANDS + [RESPONSE]
  )
  task.start()

The task starts but it fails with the following error message "Error: Image.select: Pattern 'test' did not match any bands."

When I use

featureStack = image.addBands(trainingMapReduced, [RESPONSE])

instead of

featureStack = ee.Image.cat([
  image.select(BANDS),
  trainingMapReduced.select([RESPONSE])
]).float()

I get the error message
"Error: Image.addBands: Cannot add band 'test' because it doesn't appear in srcImg.
"

The print(trainingMapReduced) looks like this:

ee.Image({
  "type": "Invocation",
  "arguments": {
    "input": {
      "type": "Invocation",
      "arguments": {
        "collection": {
          "type": "Invocation",
          "arguments": {
            "tableId": "users/xxx/xxx"
          },
          "functionName": "Collection.loadTable"
        },
        "properties": [
          "test"
        ],
        "reducer": {
          "type": "Invocation",
          "arguments": {},
          "functionName": "Reducer.max"
        }
      },
      "functionName": "Collection.reduceToImage"
    },
    "bandSelectors": [
      "test"
    ]
  },
  "functionName": "Image.select"
})

Does anybody have an idea how to solve this error?


texts: Run Python child script using Python Add-In button class?

texts:
I have a Python script which I want to execute using Python Add-In button class. Right now, I am using following:

class child(object):
    """Implementation for Toolbar_addin.button2 (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        os.system('C:/temp/child.py')

but it returns nothing when I click Python Add-In button as part of sub-menu. Is it the right way to fire child scripts using button class?

Note: child script works perfectly as a stand alone script


texts: Calculating hypsometric areas in QGIS and PostGIS

texts:
Basic Question

I am currently trying to calculate hypsometric areas of the Earth using a GeoTIFF of the ETOPO1 Global Relief map. Specifically, I require a polygon tracing all areas with a bathymetry from 0 to -X00 meters. (I require a polygon (or multipolygon) of this area for later stages in the analysis, therefore solutions for only calculating the area are inadequate).

I have attempted to do this using both QGIS and PostGIS, but I am running into a variety of problems.

My PostGIS workflow

# Import the GeoTIFF into PostGIS using raster2pgsql
raster2pgsql -s 4326 -t 50x50 -f rast -I -Y my.tif etopo_2008_raw | psql MYDB

# Generate contours of the desired area using gdal since there is no
# ST_contours type function in PostGIS that I am aware of.
gdal_contour -a ELEV -i METERS -f 'ESRI Shapefile' 'PG: dbname=MYDB host=LOCALHOST user=postgres password=PASSWORD port=PORT mode=2 schema=SCHEMA column=rast table=etopo_2008_raw' ~/PATH/etopo_contours

# Load the output shapefile back into postgres
shp2pgsql -I -s 4326 etopo_contours etopo_contours | psql MYDB

Next I use the following SQL query to 1) subset out the relvant contours, do lines to polygons then do the symmetrical difference. I apologize for how hairy this sql is.

CREATE TABLE hypsometric_area_X00 AS
WITH A AS (SELECT st_makevalid((st_dump(st_polygonize(geom))).geom) AS geom FROM etopo_contours WHERE elev='0'),
B AS (SELECT st_makevalid((st_dump(st_polygonize(geom))).geom) AS geom FROM etopo_contours WHERE elev='-X00')
SELECT ST_SymDifference(C.geom,D.geom) AS geom FROM (SELECT *, row_number() OVER (ORDER BY 1) AS rn FROM A) AS C JOIN (SELECT *,row_number() OVER (ORDER BY 1) AS RN FROM B) AS D ON C.rn=D.rn;

This query does work, but does not give a correct result. This is presumably because the st_polygonize is creating bizarre geometries from the contours. Here is a picture of the output for the 0 to -300 case. Australia is definitely not -300 meters below sea-level.

Incorrect PostGIS Output

Alternative QGIS Workflow

The QGIS workflow is conceptually identical to the PostGIS workflow.

Import the GeoTIFF using add raster layer
Use Raster->Conversion->Contour to extract the contours. (Note this uses the exact same gdal query as above.)
Subset out the contours that equal = 0 and -X00 as separate line layers.
Use Vector->Geometry Tools->Lines to Polygons to create polygons

The resulting polygons look much more sensible than what we got from PostGIS.

The last step should be relatively simple, use either Vector->Geoprocessing Tools->Difference or Vector->Geoprocessing Tools->Symmetrical Difference to find the 0 to -X00 polygon. However, this attempt consistently crashes or freezes, even if the "ignore incorrect geometries selection is chosen". I've tried it on both a Mac and Windows machine - both running QGIS 2.18.

Summary of Final Question

I am looking for any alternative workflow or modification to the current workflow to make this work in QGIS or PostGIS. (I prefer a PostGIS answer). To clarify the desired final product, I want just the blue part of picture 2 as a separate layer.

P.S.
My intern spent a week trying to get this to work in ArcGIS to no avail. I prefer a PostGIS solution because the rest of the project is coded that way, so for consistency and reproducibility it would be better.


texts: How often are MapQuest tiles updated with new OSM data?

texts: When users make edits to OpenStreetMap, how long does it generally take for those edits to show up in the MapQuest tiles?  Is this on the order of hours, days, or months?


texts: Creating vertex where two lines meet using QGIS?

texts:
I've cleaned up a road network so that I can use the shortest path tool, but some intersections don't have a vertex on one/both roads, so the result ignores this intersection.

How can I split or insert a vertex on one line where it meets another line, whether at a T-intersection or road crossing?

I've tried extending, splitting, etc. but nothing works.


texts: PostGIS (2.0?) - Native Versioning of the Datastore

texts:
I recently have been surprised to see there are not too many solutions at all discussing the best methods of versioning datastores which are made available for editing.

I am looking for a robust, efficient solution which would me to share a bunch of data with registered users and where those users can edit this geodata or associating attributes individually. My requirement though is that I need any and all edits logged so they can be rolled back to any previous state just like a wiki.

More specifically, I am utilizing a tool called CartoDB.com (built on PostGIS 2.0) through which whatever method needs to be compatible.

Very curious to learn any/all information you guys might be able to share in regards to my needs. Thanks - CH


texts: Comparing size between polygons in different UTM zones

texts:
I am going to compare size of different areas in Europe in different UTM zones. The areas I'm comparing are quite small (around 35-50 km2) and I was wondering how to do this the best way? Can I draw a polygon in one UTM zone and then move it to a different UTM zone to compare it with the polygon in the other zone? Most areas are within UTM 32N but I have some areas in UTM 29 and some in UTM 34.

The goal is to compare the size of each polygon compared to a reference zone in UTM 32.


texts: Exporting ArcGIS Online Custom Popup Window Content?

texts: Is there a way to export the content of the customized pop-ups (that uses Arcade) in a feature class to a .csv/tabular form?


texts: Polygons on Mouseover event only Activate on Borders

texts:
I have a map with several polygons that have a mouseover and click event.  The mouseover event causes the polygon to change color, and the click event brings up a custom info window.  Most of the 200 or so Polygons work, but there is a handful (about ten), that only do the mouseover or click along the borders of the polygon.

Here is the link to the map. The polygons that don't work are shown.  There is one in the middle as an example of a working polygon.

http://www.houstonportmap.com/directory.php

There are also some polylines on the map, but I've checked and the polygons appear to be correct.


texts: Cumulative sum of pixel values intercepted above a given rasters

texts:
I have 41 density calculated rasters, which raster represents a layer of a las file. I need to calculate the cumulative sum of LiDAR points intercepted above a given layer.
density_rasters a vector that includes a 45 raster layer as list. I was suggested to use Reduce function, but it creates RasterLayer with NA values. How can solve it? This is a part of for loop. If it is not clear, I can share all loops.
pcs=Reduce("+",density_rasters[(i+1):length(density_rasters)])


texts: Is it possible to draw polygons along Google Maps Polyline?

texts:
Web application has option that user draws polyline on Google Map. After he draws it, there is a button 'Add polygons', which should draw polygons along that line, based on some width and height in meters. After polygons are drawn, polyline should be removed. All polygons should have right side joined with next one (should share one side with next one). I don't know from where should I start, is it possible to draw polygon based on distance (m) from line, how to draw polygons so that they share one side?

It's clear that upper and bottom side won't be the same, that doesn't matter. Width doesn't have to be exact.

This is polyline:

myPath = new google.maps.Polyline({
                        path: flightPlanCoordinates,
                        geodesic: true,
                        strokeColor: '#FF0000',
                        strokeOpacity: 1.0,
                        strokeWeight: 4,
                        editable: true
                    });

This is fast sketch, just to get an idea what I'm talking about.


texts: Error adding MrSID raster layer to QGIS for Mac OS X 10.10.2 (Yosemite)

texts:
I want to add a MrSID layer to QGIS (version 2.6.1) with a computer running Mac OS X Yosemite. I installed QGIS through macports from the terminal and all python packages through macports, as well. I have been googling the problem and have tried to install all kinds of packages that MrSID may be dependent on, but it still doesn't work. This is the error it prompts:

  "/Users/jordi_pages/Downloads/LT_LiDAR_Sample_Data/Dallas_11SDK.sid is not a supported raster data source"

Any ideas?


texts:
I made a Python-Addin-Button with the function to create a zip. Everything works fine, but after creating the zip-file, the zipfile is locked, when the mxd-Project is open. Any suggestions why I can't open or delete the zipfile?

Here my script:

    import sys, os, arcpy, zipfile, traceback, pythonaddins
    from arcpy import *

    class CreateZIP(object):
        """Implementation for Python_Addins_addin.button (Button)"""
        def __init__(self):
            self.enabled = True
            self.checked = False
        def onClick(self):
            arcpy.env.overwriteOutput = True
            mxd = arcpy.mapping.MapDocument('current')
            mxdpath = mxd.filePath
            inWorkspace = os.path.dirname(mxdpath)
            outWorkspace = inWorkspace
            arcpy.env.workspace = inWorkspace
            outzip = zipfile.ZipFile(outWorkspace+"\Arch_Data.zip","a",zipfile.ZIP_DEFLATED)
            infolder = os.path.join(inWorkspace,"DATA")
            path = os.path.normpath(infolder)
            for (dirpath, dirnames, filenames) in os.walk(path):
                for file in filenames:
                    if not file.endswith(".lock"):
                            outzip.write(os.path.join(dirpath, file),os.path.join(os.path.basename(path),
                                                              os.path.join(dirpath, file)[len(path)+len(os.sep):]))
            outzip.close()
            z = zipfile.ZipFile("Arch_Data.zip","a")
            z.write("README.txt")
            outzip.close()


texts: White stripe when exporting QgsLayout Map to image

texts:
I have a problem when I try to export a Map in the Layout of QGIS 3 with the Python API. Here is my code :

qpt = QDomDocument()
qpt.setContent(content)
layout = QgsLayout(QgsProject().instance())
layout.loadFromTemplate(qpt, QgsReadWriteContext())
mapItem = layout.referenceMap()
layoutExporter = QgsLayoutExporter(layout)
imageSettings = layoutExporter.ImageExportSettings()
imageSettings.dpi = 150
imageSettings.generateWorldFile = 1
path = "path/to/my/iamge.tif"
layoutExporter.exportToImage(path, settings=imageSettings)

When I export my layout's map I get an image with a white stripe at the bottom like in the following file :

The width of the stripe is variable when for example a widget (like the python console) is opened or not. It seams that the map of the layout doesn't fill its whole surface because its ratio is not the same as the map canvas' one.

Is there a way to say to the layout's map to fill completely ?


texts: Convert x, y position in georeferenced image (with world file) to longitude, latitude?

texts:
What steps do I need to do in order to calculate longitude and latitude for a specific pixel in a georeferenced image? For example this image http://www.shadedreliefarchive.com/medit_Mideast_africa_copy.html contains information, .tfw file:

1627.32174969982000
0.00000000000000
0.00000000000000
-1627.32174969982000
-4957506.34276705000000
6096922.76188281000000

and .prj file:

PROJCS["Lambert Azimuthal Equal-Area",
       GEOGCS["GCS_WGS_1984",
              DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137, 298.257223563]],
              PRIMEM["Greenwich", 0],
              UNIT["Degree", 0.017453292519943295]],
       PROJECTION["Lambert_Azimuthal_Equal_Area"],
       PARAMETER["false_easting", 0],
       PARAMETER["false_northing", 0],
       PARAMETER["latitude_of_origin", 0],
       PARAMETER["central_meridian", 20],
       PARAMETER["xy_plane_rotation", 0],
       UNIT["Meter", 1]]

I've created this snippet of javascript to compute new x1 and y1 values for specific pixels in the image:

var A = 1627.32174969982000,
    D = 0.00000000000000,
    B = 0.00000000000000,
    E = -1627.32174969982000,
    C = -4957506.34276705000000,
    F = 6096922.76188281000000;

function x1(x, y) { return A*x + B*y + C; }

function y1(x, y) { return D*x + E*y + F; }

What more do I have to do to get longitude and latitude? For example C, F represents the top-left corner of the image, right? Then how do I translate C, F to get longitude and latitude?

I want do do it using javascript only, have looked at http://proj4js.org/ but don't really understand what to do with the values I get...


texts: Carto: Markers being offset or cut-off on polygon & marker map

texts:
I am seeing artifacts and markers are getting cut-off on my map. Zooming in and out can also leave an 'imprint' of what the map looked like before - either for a few seconds or until you refresh the page.

What could be causing these issues and how could I fix them?

I've added screen-grabs below.

Left: Scrambled markers/artifacts. Right: Zoom display problem

EDIT: Using marker-placement: point and marker-multi-policy: whole has almost fixed my problem. My new map is here. But some markers are still getting clipped/cut-off. How can I stop that from happening? I've tried using buffer-size as shown here but it doesn't seem to be working.


texts: Summarise rasters then project, or project then summarise?

texts:
I'm working on an analysis that involves summarising a stack of unprojected (WGS84) global 30" climate rasters (various expressions of temperature and precipitation), and generating final output in an equal-area projection.

Is it better to project the individual global rasters first (e.g. to Mollweide) and then summarise, or to summarise and then project? What factors (other than computational efficiency) should be considered in this regard?

Additionally, final output will probably be aggregated to 10 km. Is it more appropriate to aggregate raw data to, say, 0.1° and project to 10 km, or project to ~1 km and then aggregate to 10 km?


texts: Image segmentation in open-source software

texts:
For an analysis of several multilayer raster files I want to perform some kind of image segmentation (multi-resolution, ...) to extract polygon object of similar raster properties. I know that Trimble's eCognition is quite useful fot this kind of tasks. However I am trying to find a script based solution, since there are multiple (some 20s) of rasters the image segmentation schould be performed on.

I personally prefer the script languages R or Python.


texts:
I'm currently trying to create an animation with ArcScene along a given path. When I follow the Esri procedure Moving the camera or view along a path, I do not get the expected results. I end up with no data represented on my screen.

My data for this example are :

a simple DEM, with a base heights floating on a custom surface (the DEM surface)
a line with a given offset floating on the same surface. The camera is supposed to follow that line.

From this point, I select the red line and go in the Animation menu to select the "Create Flyby from Path" option. I tried to fiddle with all the options (Path source, simplification factor, Path destination, orientation settings,...). However, it seems that the camera view that does not focus on my data. It is like if the camera is moved far away from the destination data (the screen is totally white during the animation time).

When looking at the Animation manager, it seems to me that the Target:X, Target:Y and Z make sense since they do correspond to the coordinates of my data.

Yet, when I go manually to a view that I would like to obtain, the Animation Manager gives me totally different values.

Am I missing something ?

I would like to avoid the manual way to do so with the fly function as it is hard  to control precisely the flight speed/direction and impossible to reproduce this same flight several times.


texts: Having issues getting Leaflet "L.tileLayer.wms" to display Geoserver WMS data

texts:
I've set up Geoserver and WMS capabilities to host some weather data files (satellite data for example). Here's an example file:

http://xxx.xxx.xxx.xxx:8080/geoserver/opengeo/wms?service=WMS&version=1.1.0&request=GetMap&layers=opengeo:met9Qtest&styles=&bbox=-39.550275759,-33.426229830999986,60.91911924099999,67.043165169&width=512&height=512&srs=EPSG:900913&format=image/png

see code below for ip

Now, I'd like to be able to use leaflet to display weather data layers. I've been successful in getting example in Leaflet documentation to work: http://leafletjs.com/reference.html#tilelayer-wms

But when I try to do this same thing for the data on my Geoserver, the layer never shows up. Here's the snippet from my Leaflet code:

var test = L.tileLayer.wms("http://209.123.242.190:8080/geoserver/opengeo/wms",{
layers: 'opengeo:met9Qtest',
format: 'image/png',
transparent: true,
opacity: 0.7
});

I've tried changing the SRS in geoserver to EPSG:900913, 4326, 102100, 3857... still nothing. Any ideas? I can't tell if the issue I am having is on the Leaflet side or Geoserver?

Thanks,
Matt


texts:
I work with this polygon shapefile:

I can see the values 0,1 color but I can't see where the purple polygon location. In the attribute table there no more than two values, so I don't understand why the purple polygon added


texts: QGIS Categorized Symbology based on joined table

texts: I am trying to create categorized symbology for a data layer. The column that I would like to base the categories on is a result of a table join. For some reason it isn't working. All of my polygons are blue (see photo). How do I fix this?


texts: Fastest way to count the number of features in a feature class?

texts: With the introduction of the Data Access module in arcpy (30x faster search cursors), I want to know if counting features matching sql criteria is faster than the traditional MakeTableView + GetCount methodology?


texts: Choosing which attributes to return with WFS?

texts:
I am trying to use this documentation:

https://docs.geoserver.org/stable/en/user/services/wfs/reference.html#getfeature

to figure out a way where I can limit the attributes returned by a wfs call. If we use this example, I am expecting that this call would only return the GEOMETRY and the LOCALITY attribute

https://openmaps.gov.bc.ca/geo/pub/wfs?
SERVICE=WFS&
VERSION=2.0.0&
REQUEST=GetFeature&
outputFormat=application%2Fjson&
typeNames=WHSE_IMAGERY_AND_BASE_MAPS.GSR_AIRPORTS_SVW&
SRSNAME=EPSG%3A3005&
propertyName=LOCALITY&
COUNT=10


texts: How to remove a certain number of digits from the front or end of an attribute field?

texts:
Can I use the filed calculator in QGIS somehow to remove the first 4 numbers in all attribute fields in a column?

Example 16011234 -> 1234

EDIT:

Silly, I can just subtract with 16000000 in this case.

But to make it more interesting. How would I do to remove the first 2 and the 2 last numbers?


texts:
I'm in the process of defining a GML schema from a database schema, the idea is for it to be GML simple feature compliant. However I can't find a way to confirm this compliance other than reading the spec and manually checking. I'd feel happier if there was an automated tool I could pass the schema to which could check all was well. The existance of this sort of tool is hinted at in the spec (Appendix A) but I can't find one.

Does anyone know of such a tool?


texts: Edit Shapefiles Attribute data in GeoServer Programatically

texts:
I'm rendering a shapefile using GeoServer and using OpenLayers library to render it in the web browser.

Now I have created the option of editing the attribute data in web browser. I want the edited data to get reflected back in the original shapefile.

What is the best method to edit shapefile programatically using Python. I'm quite new to GIS concepts. Or do i need to convert it to PostGIS format in-order to edit the attribute data using Python?

What are the best practices and way to go?

Thanks in advance.


texts:
I'm new. I did see this was asked a couple years ago, but I didn't understand the answer. It said It's the projection. You're in ArcScene so be sure all of your layers are in Conic! I have no idea what conic means in terms or ArcGIS or ArcScene. I looked it up and found nothing.

Why do layers disappear in ArcScene when I try to drape them to the TIN?

I could set the elevation offset at 0 or 10 or 20 and it happens every time to the trees, construction and contours layers. 2 other layers remain, but they look messed up.
I saw a possible reason was the graphics card. I have an Intel UHD graphics card on this ASUS laptop.

Other specs:
Windows 10
RAM: 8 GB
Processor: Intel Core i5-8250U CPU @ 1.60 GHZ - 1.80 GHZ
64 bit OS, x64-based processor
ArcGIS 10.7


texts: Calculate Volume over sloped ground in GRASS GIS

texts:
I am trying to do some volume calculations on a DEM with GRASS GIS.

However, I believe r.volume measures everything from 0.

The ground underneath the pile I'm measuring is sloped.

I would like to calculate the volume to reflect this, but I don't know the best way to go about it.


texts: Having difficult time to implement a Rule using OpenLayers.Filter.Function

texts:
I'm trying to set up a Rule using the Function Filter.
But it does not work: the lines on that layer are just black. And more importantly, the filtering function is never executed.

Can someone please point out what is the mistake that i'm making?

Here is the code.

Thank you for your time and kind concern.

var my_filter = new OpenLayers.Filter.Function(
    function(attributes) {
       console.log(attributes);
       var x=0;
       return true;

    });
var ruleLow = new OpenLayers.Rule({
    filter:my_filter ,
  symbolizer: {pointRadius: 10, fillColor: "green",
               fillOpacity: 0.5, strokeColor: "green"}
  });

var my_style=new OpenLayers.Style( null,ruleLow);

var my_style_map=new OpenLayers.StyleMap({
    "temporary":my_style,
    "default":my_style,
    "select":my_style

});
this.vectors= new OpenLayers.Layer.Vector(
    "Vector Layer",
    {
        styleMap:my_style_map,
        renderers:this.renderer,
    }
);


texts:
I'm working with Arcgis 10.1, i'm trying to loop if condition on Python code, but they give me as error :

  Parsing error IndentationError: unexpected indent (line 18)

I don't know what is going wrong.

This is the code.

import arcpy
import math
arcpy.AddField_management("Parcelles_class_FeatureToPoi3","new_field","DOUBLE")
cur = arcpy.UpdateCursor("Parcelles_class_FeatureToPoi3")
    for row in cur:
...     PointX = row.getValue("POINT_X")
...     print PointX
...     PointY = row.getValue("POINT_Y")
...     print PointY
...     EndX = row.getValue("end_x")
...     print EndX
...     EndY = row.getValue("end_y")
...     print EndY
...     dx = EndX - PointX
...     dy = EndY - PointY
...     azimuth = math.atan2(dy, dx)
...     print azimuth
...     azimuth_gr = azimuth*200/3.14519
...     print azimuth_gr

    if (dy>=0 and dx>= 0 ):
...     azimuth_fi = azimuth_gr
...     else :
...         if (dy>=0 and dx<= 0 ):
...             azimuth_fi = 200- azimuth_gr
...             else :
...                 if (dy<=0 and dx<= 0 ):
...                     azimuth_fi = 200 + azimuth_gr
...                     else :
...                         if (dy<=0 and dx>= 0 ):
...                             azimuth_fi = 400 - azimuth_gr
...                             print azimuth_fi

    row.setValue("new_field", azimuth_fi)
    cur.updateRow(row)


texts: Formatting grid UTM co-ordinate numbers in QGIS 3.16

texts: I have used the expression @grid_number ||if(@grid_axis='x',' mE',' mN') to add the suffix "mE" and "mN" to the UTM grid labels in a layout window (e.g "550000 mE"). How do I add a comma separator for the thousands in the labels (eg convert "550000 mE" to "550,000 mE")?


texts:
I am having trouble generating a colorbar for two raster plots. I could not manage to install basemaps in my Mac so I am forced to use matplotlib.
import rasterio
from matplotlib import pyplot

april_2019rain = rasterio.open('…/rain_0419Vis.tif')

april_2020rain = rasterio.open('…/rain_0420Vis.tif')

fig, (axr, axg) = pyplot.subplots(1,2, figsize=(20,10))

show((april_2019rain), ax=axr, title='Total Precipitation April 2019').grid(True)

show((april_2020rain), ax=axg, title='Total Precipitation April 2020').grid(True)

cbar = fig.colorbar(april_2019rain)

pyplot.show()

I would really like to have the colorbar horizontally displayed at bottom of the plots.


texts: Delineating adverse/favourable areas adjacent to roads using ArcGIS Desktop?

texts:
I have a client in the forestry sector that would like me to run some maps with criteria that I’ve never managed for.
These maps will delineate (using polygons) areas of slope (0-20%,20-40%,>40%), distance from road centreline (0-125m, 125-250m, 250-350m, >350m) and “adverse” and “favourable” skids for a given “block area”.

The slope polys and distance buffers are easy enough to produce; I do this all the time.

The challenge for me will be figuring out how to identify the skid conditions. To give more insight on that, a “favourable” skid is defined as an area where harvested wood will be moved “downhill” to a road for decking. An “adverse” skid would be an area where the machine has to skid wood uphill. So this ultimately comes down to the location of the road, the +/- condition of the slope, and the aspect of the slope relative to the road.
It’s a fairly complex condition to analyze, but I’m sure it can be done. My client is asking me to replicate another GIS contractor’s version of this product… which I am fairly certain is generated with some sort of model or script.

Has anyone done anything similar, or have any insight on how I could achieve these results?

I am running ArcGIS 10.3 and use anything from 1m-20m DEM and LiDAR datasets, but each block has its own constraints.


texts: How to find the minimum perimeter oriented bounding box?

texts:
I would like to know an algorithm to calculate the minimum perimeter oriented bounding box of a point set or polygon. I know that if I use the rotating calipers, I am able to get the minimum area oriented bounding box. However, I would like the bounding box considering the perimeter, or the width.

An example of this is given in ArcGIS.


texts:
I am in the process of editing gpx files and need to write the files back to .gpx in order to put them in a GPS unit.  Basically, I am cleaning up some data and removing some track points.  However, when I write the .gpx file, it clears out the timestamp, which makes it very difficult to work with the data later.

I am using the writeOGR function from the rgdal package to export the file.  An example that can be reproduced with any .gpx file is as follows:

writeOGR(track, dsn='myfile.gpx', layer="track_points", driver="GPX", dataset_options="GPX_USE_EXTENSIONS=yes", "FORCE_GPX_TRACK=true", overwrite_layer = T)

When you read this back into R using readOGR("myfile.gpx", layer="track_points), the 'time' column is all NA.

I am looking for a way to write the original time information to the new .gpx file.

Here is an example workthrough:

library(rgdal)

track<-readOGR(system.file("vectors/test_trk2.gpx", package="rgdal"), layer="track_points")

track<-track[is.na(track$time)==F,]  #eliminate NAs from track's time column

writeOGR(track, dsn='myfile.gpx', layer="track_points", driver="GPX", dataset_options="GPX_USE_EXTENSIONS=yes", "FORCE_GPX_TRACK=true", overwrite_layer = T)

track<-readOGR("myfile.gpx", layer="track_points")  #read track back in

track@data


texts: Error DescribeData: Method areaFieldName does not exist

texts:
I am trying to Increase Polygon Area using "Increase Polygon Area" Tool in the "Production Mapping Tools" in ArcGIS 10.6, which is returning the following error: DescribeData: Method areaFieldName does not exist.

I added a new field and calculated (Calculate geometry) area and set the projection, but no luck so far.

Is there any idea how can this error be fixed?


texts: What spatial reference is used when setting Region in Grass?

texts:
I'm using Grass to process some data that is in UTM:43N. While creating the Location, it asked me for the Spatial reference, and I Entered the EPSG code for this Projection (EPSG:32643).

Now it's asking me to set the Region:

Do I enter the coordinates in EPSG:4326 (i.e. Lat Long), or in EPSG:32643 (UTM 43 N)?


texts: Seeking Random Forest R Code for Temporal Datasets

texts:
I have images from 4 different years and field datasets that correspond for each year.
Is there an R Code for Random Forest that will give me a classification map based on the 4 images and their respective field data?
I was thinking of stacking the images, but I don't think this is the right approach since each image has its own field data.


texts:
I am using GRASS 7.2.0 on Windows 8.1 and I have 107 raster files, converted from r3.to.rast which I want to display all at once.

I put in Python console (2.7.5) every name of rasters into variable d (list) and type:

for item in d:
    grass.raster.run_command("d.rast", map = item)

and it crashed if writes:

  neither monitor (managed by d.mon command) nor grass_render_immediate
  (used for direct rendering) defined

I tried to run d.mon but no help.

How can I do that?


texts: QGis Application Debugging - missing source problem at setScene method in QgsMapCanvas class

texts: I develop a app in C++ with Qgis components with VS2010 Prof on Windows 7. Compiling a release succeeded and the app´s behaviour is as expected. Nevertheless I would like to take advantage of the debug mode for further modifications. I get an an unhandled exception in the nt.dll, missing source, whenever the debugger comes to the setScene(mScene) method. I made a folder with all available '*.cpp' files from the QT4 and the QGis source files, so all sources should be available. Does anybody know which source file might be missing at this point?


texts: Merging a temperature point layer with a vector layer for statistical analysis in QGIS

texts:
I have a point layer with x,y,z coordinates with geothermal temperature data. This data should be related to a given vector layer of political communities. In order to analyse every single community to their temperature distribution according to the depth (min, max, mean, standard deviation).

I tried it with the Spatial data processing framework plugin "merge attributes by position" (not sure if its the right translation from German). With this I get the sum, min, max results for each community.

I would like to get all the temperature data within every single community, to have further possibilities of analysing the data.


texts: How to use tiles in this strange tiling scheme?

texts:
I'm looking at cell coverage maps such as this one, and trying to figure out how to get a GIS program to use it as a source of information. Ultimately, I need to be able to combine the output with other cell coverage maps from different sources. But in order to do that, I need to understand how each one works and how to align them with each other.

Strangely, each of the five zoom levels offered has a different tile size:

160x160 pixels
128x128
145x145
198x198
182x182

I've never seen a map tile service like that before. Usually the tile sizes are all the same. Do map-service-consuming tools (I'm thinking of open-source tools like Tilecloud) have the ability to consume maps whose layers have different - practically random - tile sizes?

Additionally, I'd like to know how the tiles are georeferenced. Maybe that's not exactly the right way to say it ... but I need to be able to combine tiles from different maps so that they're aligned properly.

Maybe I could get the information needed for that if I knew what kind of map tile API is being implemented by http://gsma.streetmap.co.uk/custom/barts/gsma/. Maybe they have a REST URL that gives data about spatial reference and so on. But I haven't been able to figure out what API they're implementing.

Using the browser debugging tools, I see URLs like http://gsma.streetmap.co.uk/custom/barts/gsma/gmap.srf?o=1&i=502/951/Zoom3/Zoom3_502_Tile007_9511638.png
for the tiles.

The mobileworldlive.com page loads several javascript files, including
two from gsma.streetmap.co.uk:
http://gsma.streetmap.co.uk/custom/barts/gsma/SMJ_GSMA.srf?id=502&layer=951&nd=http://mobileworldlive.com/maps/images/map_na.jpg
and
http://gsma.streetmap.co.uk/SMJ_ObsEmo.js
(The first one has an "include" that loads the second.)

They contain a bunch of custom map display code. I suppose I could try to learn about the API through reading that code, but that seems like a time-consuming and error-prone way to do things if this is implementing an API that's documented somewhere.

Thanks for any suggestions on how to attack this.


texts: Converting lat,lon to x,y coordinates with pyproj

texts:
I am trying to convert lat,lon to x,y coordinates in France, using the pyproj library and what I am doing is obviously wrong.

Here is some of my code:

import pyproj
import math

PROJ='+proj=utm +zone=31, +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

def LatLon_To_XY(Lat,Lon):
    p1=pyproj.Proj(PROJ,preserve_units=True)
    (x,y)=p1(Lat,Lon)
    return(x,y)

def XY_To_LatLon(x,y):
    p1=pyproj.Proj(PROJ,preserve_units=True)
    (lat,lon)=p1(x,y,inverse=True)
    return(lat,lon)

def distance(x1,y1,x2,y2):
    d1=x1-x2
    d2=y1-y2
    out=math.sqrt(d1*d1+d2*d2)
    return(out)

Problem is when I am validating this computing the distance between Paris and Toulouse, I get 816km where the distance is supposed to be around 580km.

I guess I have wrong parameters for the projection definition.

Here are the values displayed:

For Paris lat,lon (48.856614, 2.3522219) and x,y: 6253348.594037977 374425.11288885673
For Toulouse lat,lon (43.604652, 1.444209) and x,y: 5453729.038474649 210744.18988320036


texts: Calculating minimum distance using "Distance Matrix" in QGIS

texts:
I want to calculate minimum distance between points and a line in QGIS. I use "Distance Matrix", point shapefile contains 1337 points data, and I made unique ID as sequence number. A line is actually coast line and I convert text file to vector, then extract nodes using "Extract Nodes".

I followed manual and tutorial, but the output is different with what I expected. They mentioned that if I want to calculate minimum distance note check k=1, but I think my result contains the farthest distance. The output file's target ID is all the same. And I wonder what is the exact meaning of unique ID?

Please check and advise about my actions.


texts: ArcMap Random Points Minimum Distance Field not able to use units

texts:
I would like to create random points with a minimum distance based on an attribute field (polygons with a larger area need random points with more space between them).

In theory this can be done. The random points tool has such a feature. But it doesn't recognize a text attribute. The buffer tool has a similar feature, and recognizes if you put the number + "Miles" in a text attribute column.

Random Point minimum distance only recognizes number attribute columns, and it is not clear which units it is converting those distances to. I imagine they should be in meters, for my data and data frame are both in British National Grid, yet points are created in polygons closer than what is stated in the field.

Perhaps I'm missing something.


texts: PRACH data plotted in a better way but not completely

texts:
I am close to finishing the visualization of PRACH in MapInfo.

Here's my output so far (for 4 counters)

Unfortunately, for some sectors it is not working as it should be.

Can you please help out?

Here's my code:

Select * from PRACH_data where Obj into ToUpdate

Fetch First From ToUpdate
Do while not EOT(ToUpdate)
   dim x as float, y as float
   azimuth = ToUpdate.Azimuth
   BW = ToUpdate.HBW
   s = azimuth -(BW/2)
   e = azimuth + (BW/2)
   latitude = ToUpdate.Latitude
   longitude = ToUpdate.Longitude

   circ = CreateCircle(ToUpdate.Longitude, ToUpdate.Latitude, ToUpdate.Dist_Class3)
   Create Arc into Variable pie
      (ObjectGeography(circ, 1), ObjectGeography(circ, 2))
      (ObjectGeography(circ, 3), ObjectGeography(circ, 4))
      s e
   pie=converttopline(pie)
   alter Object pie
      Node Add (ToUpdate.Longitude, ToUpdate.Latitude)
   pie=converttoregion(pie)
   Insert into WindowInfo(FrontWindow(),10) (obj) Values (pie)

   circ = CreateCircle(ToUpdate.Longitude, ToUpdate.Latitude, ToUpdate.Dist_Class2)
   Create Arc into Variable pie
      (ObjectGeography(circ, 1), ObjectGeography(circ, 2))
      (ObjectGeography(circ, 3), ObjectGeography(circ, 4))
      s e
   pie=converttopline(pie)
   alter Object pie
      Node Add (ToUpdate.Longitude, ToUpdate.Latitude)
   pie=converttoregion(pie)
   Insert into WindowInfo(FrontWindow(),10) (obj) Values (pie)

   circ = CreateCircle(ToUpdate.Longitude, ToUpdate.Latitude, ToUpdate.Dist_Class1)
   Create Arc into Variable pie
      (ObjectGeography(circ, 1), ObjectGeography(circ, 2))
      (ObjectGeography(circ, 3), ObjectGeography(circ, 4))
      s e
   pie=converttopline(pie)
   alter Object pie
      Node Add (ToUpdate.Longitude, ToUpdate.Latitude)
   pie=converttoregion(pie)
   Insert into WindowInfo(FrontWindow(),10) (obj) Values (pie)

   circ = CreateCircle(ToUpdate.Longitude, ToUpdate.Latitude, ToUpdate.Dist_Class0)
   Create Arc into Variable pie
      (ObjectGeography(circ, 1), ObjectGeography(circ, 2))
      (ObjectGeography(circ, 3), ObjectGeography(circ, 4))
      s e
   pie=converttopline(pie)
   alter Object pie
      Node Add (ToUpdate.Longitude, ToUpdate.Latitude)
   pie=converttoregion(pie)
   Insert into WindowInfo(FrontWindow(),10) (obj) Values (pie)

   Fetch Next From ToUpdate
Loop

Commit Table ToUpdate Interactive

Also, how can I calculate the area with every radius I'm plotting the arcs? in the table, I don't have "Area" as a column but when I click on a polygon, it is giving me a value for the area.. Maybe, if I have the Area column, I can order them in a descending order and then map them... that should work, what do you guys think?


texts: Convert spatial table into non-spatial table in Oracle

texts:
I want to convert a spatial table, registered with Geodatabase (SDE), into a non-spatial table using SQL Developer, what would it be the most efficient/elegant way of doing it?

Obviously, first I would make a copy of it before changing anything.

Having said that, I would also like to maintain the foreign keys, indexes and constraints of the non-spatial columns.

which would it be the best option from below?

DELETE the 'SHAPE' and the 'SE_ANNO_CAD_DATA' columns.
DROP the Spatial Metadata/Spatial indexes.
Export the table into a text/csv/excel/etc file and upload it as a new table without spatial attributes and then, re-create the foreign keys, etc.


texts: How do versioned views work with multiple editors?

texts:
I understand that in the database environment, by executing

version_util.set_current_version('version1')

Editor1 is able to write SQL statements against [table_name]_vw views to edit the data at version1.

But if Editor2 (from another workstation) executes

version_util.set_current_version('version2')

the [table_name]_vw views should be populated with the data at version2, right? In this case, if Editor1 run a query against the [table_name]_vw, will he still get the data at version1? Or is he still editing the right version?


texts: Only label visible attributes using ArcGIS Pro

texts:
I have a huge polygon dataset of land boundaries of which I only have 40 land boundaries visible.
Is there a way to only label these visible polygons in ArcGIS Pro?


texts: Cannot transform point in OpenLayer (not able to find proj4js)

texts:
I am no more able to reproject a point in OpenLayers 6:
let center = ol.proj.transform(
  [2600000.0, 1200000.0],
  'EPSG:2056',
  'EPSG:4326'
);

=>

But this is working:
let center = ol.proj.transform(
  [46.88668, 7.47806],
  'EPSG:4326',
  'EPSG:3857'
);

Adding this line doesn't work: ol.proj.proj4.register(proj4); (as suggested in this GitHub page: https://github.com/openlayers/openlayers/issues/9199):
Even when defining a proj4js string before, nothing changes:
proj4js.defs["EPSG:2056"] = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs";
ol.proj.proj4.register(proj4);

Is OpenLayers not able to natively handle other standard projection than EPSG:4326 and EPSG:3857?
What did I missed?
OL source that I use is from here: https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.5.0/build/ol.js
And I proj4js from here: https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.7.0/proj4.js


texts: Changing attribute values in a memory layer

texts:
I'm trying to change attribute values in a memory layer. My memory layer is defined in this way:

pro.addAttributes([QgsField("nom", QVariant.String), QgsField("valeur", QVariant.Int)])

I made several tries:

attrs = { 0 : "hello"}
bFlag3=vl.dataProvider().changeAttributeValues({ fet.id() : attrs })

Or in this way :

bFlag3=vl.changeAttributeValue(fet.id(),0,'hello')

Each time my boolean flag returns True but I never see any change in my memory table. I make a

vl.commitChanges()

in both cases.

When I try these two pieces  of code on the two other layers (which aren't memory layers...) I have in my project, I see the attributes change...

Can anyone help me troubleshoot this?


texts: stack 10 modis ndvi geotiffs and scale them by 0.0001

texts:
I have 10 geotiffs Modis ndvi named 2005.tif ... 2014.tif
I want to stack them all and create timeseries. i know how to do that, but the stacked images still need scaling by 0.0001. how can i bulk scale them using gdal -scale function before stacking them?

Tried the first line on OSGeo, but when openning the new tifs all I can see is 0s and 1s. any idea?


texts: Enabling Jupyter nbextension for ipyleaflet?

texts:
I was trying to test out this script https://github.com/tylere/PyDataNYC2017/blob/master/ipynb/satellite_analysis.ipynb but got stuck at the iteractive mapping part.

import ipyleaflet
map1 = ipyleaflet.Map(zoom=3, layout={'height':'400px'})
dc = ipyleaflet.DrawControl()
map1.add_control(dc)
map1

It produces something like this:

Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …

ipywidgets package also produces something similar:

VBox(children=(Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attr…

Then I tested it using jupyter nbextension list but it shows me:

Known nbextensions:.

Then I try to fix it according to this post https://github.com/jupyter-widgets/ipyleaflet/issues/173

jupyter labextension install @jupyter-widgets/jupyterlab-manager
jupyter labextension install jupyter-leaflet

However, in Command line it shows me

Error executing Jupyter command 'labextension': [Errno 2] No such file or directory

I also tried this solution Ipyleaflet map object doesn't display in Jupyter Notebook but it gets created

jupyter nbextension enable --py --sys-prefix ipyleaflet

shows me error messages:

Traceback (most recent call last):
  File "/home/gjiang/.local/bin/jupyter-nbextension", line 11, in <module>
    sys.exit(main())
  File "/home/gjiang/.local/lib/python3.6/site-packages/jupyter_core/application.py", line 266, in launch_instance
    return super(JupyterApp, cls).launch_instance(argv=argv, **kwargs)
  File "/home/gjiang/.local/lib/python3.6/site-packages/traitlets/config/application.py", line 658, in launch_instance
    app.start()
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 988, in start
    super(NBExtensionApp, self).start()
  File "/home/gjiang/.local/lib/python3.6/site-packages/jupyter_core/application.py", line 255, in start
    self.subapp.start()
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 896, in start
    self.toggle_nbextension_python(self.extra_args[0])
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 872, in toggle_nbextension_python
    logger=self.log)
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 483, in enable_nbextension_python
    logger=logger)
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 386, in _set_nbextension_state_python
    for nbext in nbexts]
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 386, in <listcomp>
    for nbext in nbexts]
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/nbextensions.py", line 351, in _set_nbextension_state
    cm.update(section, {"load_extensions": {require: state}})
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/config_manager.py", line 136, in update
    self.set(section_name, data)
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/config_manager.py", line 110, in set
    self.ensure_config_dir_exists()
  File "/home/gjiang/.local/lib/python3.6/site-packages/notebook/config_manager.py", line 67, in ensure_config_dir_exists
    os.makedirs(self.config_dir, 0o755)
  File "/usr/lib/python3.6/os.py", line 210, in makedirs
    makedirs(head, mode, exist_ok)
  File "/usr/lib/python3.6/os.py", line 210, in makedirs
    makedirs(head, mode, exist_ok)
  File "/usr/lib/python3.6/os.py", line 220, in makedirs
    mkdir(name, mode)
PermissionError: [Errno 13] Permission denied: '/usr/etc'

I'm using python 3.6.7 (default, Oct 22 2018, 11:32:17) on jupyter notebook on Ubuntu 18.04


texts: Using multiple Calculate Value tools in Model Builder

texts:
I am trying to run a model on a dataset, and in doing so I create new features from the input spreadsheet. The NAME_SHIP field is being used to create the names of the new feature classes but the field is set up with illegal characters. I am already using one CALCULATE VALUE tool in my model but found that some of the ship names also include an " ' " (apostrophe) in them. Is there a way to either:

Include an expression in the already existing CALCULATE VALUE to rid the names of the apostrophe? or

Where do I insert another CALCULATE VALUE tool in the model to fix for both issues while creating a new FC?


texts: How to convert a shapefile into OSM .PBF format

texts:
I have a large .gdbtable file (1.3 GB) with a road network that I would like to convert to .PBF format keeping the topology information.

I can load the .gdbtable file in QGIS so I thought it should be straightforward to write/export the file in .PBF but I cannot find this option. I would rather use QGIS to do this but an ArcGIS solution would be welcomed as well.

Any ideas on how to do this?


texts: Appending feature classes from one geodatabase to another geodatabase for python

texts:
I am attempting to write a stand alone python script which will loop through the first GDB and append each feature class to its corresponding feature class in the second GDB.  The error code that I keep getting is:

  "ExecuteError: Failed to execute. Parameters are not valid.  ERROR
  000338: Inputs must be either all Feature Classes, Tables or Rasters;
  not mixed.  Failed to execute (Append)."

I believe this has to do with reading the pathways to the GDB.  It is important to note that with in the second GDB the target feature class is within a feature dataset.  My arc version is 10.2 if this matters.

import arcpy
import os
import traceback

#Set up workspace
outworkspace = "C:\\Users\\bjbart\\Desktop\\Marathon_Loading\\SURVEY_POINTS_V10_1.gdb"
workspace = "C:\\Users\\bjbart\\Desktop\\Marathon_Loading\\out.gdb"

arcpy.env.workspace = workspace

#loops through the first GDB
featureclasses = arcpy.ListFeatureClasses()

for fc in featureclasses:

    shp_file = fc
    print fc

    if shp_file == "BEND":
        arcpy.Append_management(os.path.join(workspace, "shp_file"), "C:\\Users\\bjbart\\Desktop\\Marathon_Loading\\SURVEY_POINTS_V10_1.gdb\\PIPE_FEATURE\\BEND", "NO_TEST", "","")

        print('BEND File Has Been Appended To Target Database')

    else: print('Error shape file did not get read')'


texts: Handling OpenStreetMap 64-bit identifiers correctly in QGIS?

texts:
An OpenStreetMap wiki states, "QGIS will no longer be usable with newer OSM data spreading all over the world" because it's "incompatible with OSM data containing 64-bit Identifiers."[1]  The wiki claims "QGIS may still be used for retrieving data however it is possible for the geometries to look very incorrect, with ways mixing old nodes below the 32-bit limit and new nodes not found after they are truncated to 32-bit or at wrong positions: data retrieval may unpredicatably give inconsistent results at any time."[2]

Does QGIS support 64 bit identifiers?

Has this been updated with the newer QGIS version, is it currently being updated, or is there some work-around to avoid issues with newer data?

[1] http://wiki.openstreetmap.org/wiki/64-bit_Identifiers
[2] http://wiki.openstreetmap.org/wiki/QGIS


texts: Querying with Large FeatureCollection in Google Earth Engine using Python?

texts:
I am trying to process a Tanzania shapefile downloaded from GADM Global Administrative Areas.

    # im -> {Image} ee.Image({...})
    # self.geom_coll -> {FeatureCollection} ee.FeatureCollection({...}). containing
    # 3000 features.
    # spacereducer() -> ee.Reducer.mean
    # self.scale -> 10 #Changing this value to small number gives error

    feats = im.reduceRegions(self.geom_coll, spacereducer(), self.scale)
    flist = getInfo_werrorcontrol(feats,
                          self.errorcheck)['features']

.

def getInfo_werrorcontrol(featureCollection, errorcontrolon=True):
    """
    Wrapper to add error control to GEE evaluations.

    For large computations GEE sometimes times out and needs to be
    restarted. This does so in a controlled manner with out
    interrrupting the program flow.
    """
    if errorcontrolon:
        i=0
        while True:
            try:
                with timeout.timeout(10*60):
                    return featureCollection.getInfo() # In this line I am getting exception.
            except NameError:
                exc_type, exc_value, exc_traceback = sys.exc_info()
                lines = traceback.format_exception(exc_type, exc_value, exc_traceback)
                print ''.join('!! ' + line for line in lines)
                i+=1
                print 'attempts: '+str(i)
                if i > 20:
                    raise ValueError('to many attempts')
                elif i > 10:
                    print 'waiting 2 minutes'
                    time.sleep(60*2)
    else:
        return featureCollection.getInfo()

Exception:-
Changing self.scale to 10 gives me following error for line: featureCollection.getInfo()

  ee.ee_exception.EEException: Server returned HTTP code: 413

Changing self.scale to 1000 gives me following error:

  ee.ee_exception.EEException: Computation timed out

What is the correct way to process shapefile with the larger region?


texts: Calculating distance from geolocation in Java

texts:
I'm trying to measure distance from one geolocation to another in Java. On stackoverflow I came across the following post and google developer page, giving me the following formula.

6371* acos(
cos( radians(37) )
* cos( radians( lat ) )
* cos( radians( lng ) - radians(-122) )
+ sin( radians(37) )
* sin( radians( lat )
) )

I created the following Java code.

private static double distanceCalc(float lat, float lng) {
    return 6371 * Math.acos(
            Math.cos(40.84916900104506 * Math.PI / 180)
                    * Math.cos((lat * Math.PI / 180))
                    * Math.cos((lng * Math.PI / 180) - (-73.86835600032798 * Math.PI / 180))
                    + Math.sin((40.84916900104506 * Math.PI / 180))
                    * Math.sin((lat * Math.PI / 180)));
}

but it's giving me results that are way off.  The 40.84916900104506, -73.86835600032798 geolocation is in New York but I get a value returned of around ~6090 kilometers. First I thought I just entered the lat / lng in the wrong variables but that doesnt seem to be the case.

I couldnt use Math.toRadians as it asks for an integer and not a float, which will be too inaccurate so lat * Math.PI / 180 should do right? Math isn't exactly my strongest field.