texts: Leaflet: 4 overlay layers, 1 click, get properties of all intersecting polygons.

texts:
I have a leaflet application where I have 4 layers that have many intersecting polygons such as, a counties layer, a school districts layer,  a cities layer, and a layer of miscellaneous areas. All cities and miscellaneous areas physically intersect with counties and school districts. Some miscellaneous areas intersect with cities.  I need "on" my click to get a popup that tells me the feature properties of each polygon at each layer - if that layer is indeed selected.  I am having a terrible time working on a solution.

I love how I can loop over each layer in the leaflet L.featureGroup function.  But I am stumped with querying each layer's features past that. I see the ._layer elements returned in the console.  I suspect I need to detect intersections. Can someone point me in the right direction in terms of plugins and or leaflet methods.


texts:
The objective is to create a global geotiff raster at 30s resolution (appr 1km) with the pixel value the area of that pixel. I have multiple tools but struggle to find the easiest way. I am aware of this question  and probably prefer an answer using GDAL.

QGIS:

create grid,
use raster calculator to obtain area per polygon
rasterize featurr\es

Matlab:

Areamat function

Earth Engine:

ee.Image.pixelArea()

ArcMap:
?

GDAL/ Python
?

Other
?


texts: Dealing with empty values in xml using FME XML Templater?

texts:
I am converting excel to xml using FME 2016 and the XML Templater.

When a field is empty, my xml is invalid because my schema requires

<TAGNAME xsi:nil="true" />

but instead gets

<TAGNAME />

which doesn't validate against empty values.   Does anyone know how FME can resolve this?


texts:
Is there a way to change the modifier key for the mouse wheel zoom in smaller increments?

I'd like to change it from Ctrl to Shift because I mostly use other software that works this way so I always hit the wrong modifier in QGIS and can't get used to it.

I looked in the Keyboard Shortcuts, but didn't find it. This menu seems to be only for keyboard shortcuts, not mouse interaction.


texts: Use ArcPy/Python to print out the number of selected features based on different selection criteria

texts:
I'm working with Parcel feature class in ArcGIS Pro 2.2 and trying to figure out the number of parcels that fit into each category. And each category consists of two criteria:

The development ratio of the total parcel
The acreage of the parcel

The input data table:

The circled fields are the criteria I am using. And the specific classes are shown in the result spreadsheet screenshot below.

This is the result I generated:

I ran the arcpy script for each individual category:

arcpy.SelectLayerByAttribute_management('ParcelBoundaries_DevTest', 'NEW_SELECTION', '"Developed" > 0.25 and "Developed" <= 0.5 and "GISACRES" <= 1')

arcpy.SelectLayerByAttribute_management('ParcelBoundaries_DevTest', 'NEW_SELECTION', '"Developed" > 0.5 and "Developed" <= 0.75 and "GISACRES" <= 1')

arcpy.SelectLayerByAttribute_management('ParcelBoundaries_DevTest', 'NEW_SELECTION', '"Developed" > 0.75 and "Developed" <= 1 and "GISACRES" <= 1')

After each run, I will see the number of the selected features in the attribute:

Then I will manually input the number into the spread sheet. This is still very tedious. Is there any way that I can use Python or ArcPy to generate the number of the selected features?


texts: Exporting shapefile to PostGIS database in QGIS

texts: How do I export a vector shapefile to PostGIS database? I am asking what format I have to export it.


texts:
I have a polyline spatialite database, tried to perform repair geometry in QGIS using python with the following code:

from pyspatialite import dbapi2 as db

conn = db.connect(r'D:\db.sqlite')
cur = conn.cursor()

repairGeometry = "UPDATE the_table SET GEOMETRY = ST_MakeValid(GEOMETRY) WHERE ST_IsValid(GEOMETRY) = 0;"
rs = cur.execute(repairGeometry)

I am getting this error:

rs = cur.execute(repairGeometry)
OperationalError: no such function: ST_MakeValid


texts:
There are numerous commands and libraries that I see referenced or prefixed with OGR. What does the OGR stand for?

ogrinfo
ogr2ogr


texts: Methods for determining landing hazard for NASAs Perseverance Rover

texts:
I was wondering if anyone knows the methods that were used to create the bottom map in the figure below, which is to indicate landing hazard for the Perseverance Rover. Slope and ruggedness are a couple of variables that might have been employed, but it would be nice to have a journal article or technical report that goes into more detail.
Here are links to where the map is used and landing hazard is briefly discussed, but no further sources seem to be provided:
https://www.nasa.gov/image-feature/jpl/avoiding-hazards-at-jezero-crater
https://www.jpl.nasa.gov/images/jezeros-hazard-map


texts: Finding datasets embedded in R GISTools package?

texts:
How can I find the embedded datasets in R package GISTools?

I can find examples of datasets (Georgia, Newhaven etc.), but I cannot see exactly which datasets are available.


texts: Numpy error when computing scores based on percentiles

texts:
I  want to calculate scores in one field based on the percentile (in steps of 10) the row falls into in another field. I found a solution here (Calculating Percentiles in ArcMap?) that I tried out, but I now get the error:

unsupported operand type(s) for *: 'numpy.ndarray' and 'numpy.ndarray'

My code looks as follows:

input = "Buffer"
arr = arcpy.da.FeatureClassToNumPyArray(input, ('Shape_Area'))

p1 = np.percentile(arr, 10)
p2 = np.percentile(arr, 20)
p3 = np.percentile(arr, 30)
p4 = np.percentile(arr, 40)
p5 = np.percentile(arr, 50)
p6 = np.percentile(arr, 60)
p7 = np.percentile(arr, 70)
p8 = np.percentile(arr, 80)
p9 = np.percentile(arr, 90)
p10 = np.percentile(arr, 100)
#use cursor to update the new rank field
with arcpy.da.UpdateCursor(input , ['Shape_Area','ConS']) as cursor:
    for row in cursor:
        if row[0] < p1:
            row[1] = 0.1  #rank 0
        elif p1 <= row[0] and row[0] < p2:
             row[1] = 0.2
        elif p2 <= row[0] and row[0] < p3:
             row[1] = 0.3
        elif p3 <= row[0] and row[0] < p4:
             row[1] = 0.4
        elif p4 <= row[0] and row[0] < p5:
             row[1] = 0.5
        elif p5 <= row[0] and row[0] < p6:
             row[1] = 0.6
        elif p6 <= row[0] and row[0] < p7:
             row[1] = 0.7
        elif p7 <= row[0] and row[0] < p8:
             row[1] = 0.8
        elif p9 <= row[0] and row[0] < p10:
             row[1] = 0.9
        else:
             row[1] = 1

        cursor.updateRow(row)

When searching on the internet I found that the problem behind this error often lies in the values being strings rather than floats or integers so I edited my code to read:

p1 = np.percentile(arr, 10, np.float)
p2 = np.percentile(arr, 20, np.float)
p3 = np.percentile(arr, 30, np.float)

And now I got the error message:

'type' object is not iterable

What can I do to solve this issue, or is there maybe another way to calculate the scores?


texts: Problems converting GTiff to PDS with gdal

texts:
I'm trying to convert a GTiff raster to a PDS raster (NASA Planetary Data System) using gdal_translate with this command:

gdal_translate -of PDS -ot Float32 input.tif output.IMG

But I get this error message: 0ERROR 6: GDALDriver::Create() ... no create method implemented for this format.

Did I use the wrong input or is GDAL Drivers missing? I have GDAL 2.1.0


texts: Finding least cost path to diffuse resource using ArcGIS Pro?

texts:
I am trying to solve a variant of the least cost path analysis problem (http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/creating-the-least-cost-path.htm) where instead of finding the shortest path between two points based on some cost surface, I want to find the area from which a resource can be extracted based on some cost surface. e.g. imagine a few cities in a region with demand for wood from nearby forests. each city will cut down trees starting with those closest to them.

Is there a way to model this in a GIS context?


texts:
I'm generating an image from QgsComposer composition through QgsLayoutExporter.exportToImage() method that allows saving the image to disk.

I wonder if there is any way to get image binary data into a variable?

I need to send this image later to a server to process, so I don't need a physical file at all. One way would be saving the file in /tmp and then deleting, but it requires unnecessary interaction with HDD and looks ugly.

I tried to find another solution using eg. generating an image from mapCanvas, but similarly, they provide only saving data to a file.

I'm using  QGIS API ( 3.3.0-Master)


texts:
I am trying to import the map on top of the page - into QGIS
https://www.licadho-cambodia.org/land_concessions/
I am reading the page source but not sure what is the URL I need.
I would like to import the map and the labels in the map.


texts: Migrating code to 10.1 and recycled search

texts:
I am migrating 9.3 Desktop extension to 10.1

I am getting this bug when I am trying to update a row in the geodatabase

"Cannot call Store on a recycled row while editing"

And seems the reason is, I am doing recyclable search like this

eatCursor = featLayer.Search(spatialFilter, true);

So, I changed it to false recyclable.
And it worked.

But my question is why 10.1 is not able to update recycled row, where 9.3 was able to?
what the reason of this change?
Is there any reason 10.1 disabled this feature?


texts: Find connected features by type in utility network?

texts: I have a storm drain geodatabase with digitized direction flow.  I am trying to find catch basins (upstream) that drain to other catchbasins (downstream) as opposed to a junction structure or an outlet.  The purpose is to be able to add trash collection devices on downstream catch basins instead of every catch basin.  So I need to find where these CBs are connected to each other.  Is their a programmatic way to do this?


texts: Auto-populate field based on other field when selected

texts:
I’m dealing with thousands of points in my project and requires me to select points base on distance from other line features (i.e roads). Since the possibility for some points to overlap are high, I design my work so that the selected points will turn red whilst non-selected points remains black. I don't want the points to be selected twice and belongs to two roads.
For this to work, I added another 2 text field i.e “Road” and “Marked”. After selecting the points base on location from the road, l then select all the points and fill the “Road” field with new road name, and on “Marked” field as X to change its color and differ it from non-selected points. By this, I can later remove the selected points using ‘base on attribute’ option.
There will be hundreds of road to deal with and be given new name but no matter what its name, the “Marked” field will be an X and the points that corresponding to it will change color. I can manually fill the the “Marked” field using the field calculator but am looking for the “Marked” field to be filled automatically every time the “Road” field is not empty, or else it will remain blank.
I'm using ArcMap 10.6.


texts:
I have a python script in the field calculator that reads dates and populates a status field. When there is an empty value in the "actual" field, that is the status of the project. Here is the code that I use:

def cal(rec_a, soc_a, stk_a, plat_a, perm_a, per_a, status):

    if rec_a == " ":
      return "Field"

    elif soc_a == " ":
      return "Consult"

    elif stk_a == " ":
      return "Stake"

    elif plat_a == " ":
      return "Plat"

    elif perm_a == " ":
      return "Permit"

    elif per_a == " ":
      return "Waiting"

    elif (per_a is not None) or (per_a == " "):
      return "Permit Rec"

And the expression is:

cal (!Recon_act!, !SOcons_act!, !Stake_act!, !Plat_act!, !permit_act!, !PerRec_act!, !status!)

How do I use that same expression, but instead of reading empty fields it reads Null values?

if rec_a IsNull: ?

if rec_a == Null: ?


texts: Measure shortest distance around polygon

texts: How do I measure the shortest distance between two points (1 and 2) with respect to a polygon? So the shortest way around the polygon from 1 to 2.


texts: Clip a NetCDF file using a shapefile with Python

texts:
I have global net radiation data but I want to extract for my country.  I've loaded the NetCDF but I don't know how to proceed further. The following is a similar question
(Extract data from netcdf using a shapefile polygon in python)  but the answer talks about running from command line but I'm using Spyder.

If someone could help me with the code, I know that I need to use some command like clip or extract but I'm not sure how to.


texts: How to Create a Version 10 File Geodatabase from an XML from ArcGIS 9.2?

texts:
A client has given me an XML document for the schema of an ArcGIS 9.2 Geodatabase.  The final goal will be to create a version 10 version of the geodatabase but the client doesn't have version 10 yet. I have been trying to create a new geodatabase using this XML file but am not having any luck.

So far, I have tried to create a new 9.2 Geodatabase using the Create File GDB Tool.  Once I have created the version 9.2 geodatabase, I tried importing the XML file by right clicking on it in ArcCatalog.  Every time I import it, I get the following error:

I've tried importing data as well as schema only.  I've also tried importing this xml file into a version 9.3 and 10.0 geodatabase.  All with the same result.

XML file is about 7MB.

Does anyone know why this might be happening?

EDIT:

I tried running the Import XML Workspace Document tool and it displays this error.  Is the XML file corrupt?  How could this have happened?

EDIT2:

The client has sent me a different XML file (the first one was opened in notepad and she thinks she may have corrupted it).  However, I'm getting a seperate error now.  I'm getting the following:

I Googled the above error and it may have something to do with the fact that the client is using ArcFM.

In answer to the comments, I've opened the xml file in Internet Explorer and Notepad++ and it is definitely an XML Workspace file.  It has over 100,000 lines so I wasn't keen to read it all, but wehre it says version it does say "sde.DEFAULT".  I don't have SDE set up in my office, could that be the issue?

Should the client export the data first to a file geodatabase?


texts: Add rasterlayer with defined symbology in ArcPy

texts:
I am trying to add a raster layer in ArcPy with arcpy.mapping.Layer(Raster). By default, the layer is added with "Stretched" symbology. However, for further use of the raster, I need "Unique Values" instead.

Is there a way to directly add a Layer with the predefined symbology  "unique value" to the default map document in ArcPy?

So far I have tried to change the symbology with symbologyType, but it does not seem to work.

arcpy.gp.Reclassify_sa(ReclassifyFilteredAccumulation, "Value", "1 1;NODATA 0", ReclassifyVectorize, "DATA")
VectorizeLayer = arcpy.mapping.Layer(ReclassifyVectorize)
if VectorizeLayer.symbologyType == "UNIQUE_VALUES":
    VectorizeLayer.symbology.valueField = "VALUE"
    VectorizeLayer.symbology.addAllValues()
arcpy.mapping.AddLayer(df, VectorizeLayer,"BOTTOM")


texts:
I have a dataset in the from:

gps.accuracy    latitude    longitude   timestamp   user
37.959  -9.6580833  -35.7001122 1434057621378   5573502c150000c10136e51b
37.959  -9.6580833  -35.7001122 1434059660236   5573502c150000c10136e51b
37.959  -9.6580833  -35.7001122 1434060748144   5573502c150000c10136e51b
37.959  -9.6580833  -35.7001122 1434057621378   5573502c1000000000000000
37.959  -9.6580833  -35.7001122 1434059660236   5573502c1000000000000000
37.959  -9.6580833  -35.7001122 1434060748144   5573502c0000000000000000

There are two user (this is the id and several points for each one of them).

IN QGIS 3 I manage to add the points on map using the Data Source Management, Delimited Text. This is an example (of course above dataset is a sample):

I wonder if it is possible for every user, to have a different color in order to see the different points for the different user. Also is it possible to draw a line between points of the same user?


texts: Is it possible to export a WMS-layer to dwg or dxf. I would like to use lines, circles etc. in Autocad, but I can't find out if it is possible?


texts: Looking for "big names" that use open source WMS (like geoserver / ms4w / qgis server) for presentation

texts: I am looking for famous companies / "big names" (like army , goverments , city halls) that chose to use open source map servers and not paid solutions like Esri and such for a presentation I am doing.


texts:
Can we change the view of the map that comes from MapStore in GeoNode?
How to customize the panel that shows the result of identify action to be a popup instead of the right-side panel?
I have tried to check the map app in GeoNode but still couldn't resolve it.
I couldn't understand how {% get_map_new  %} works.


texts: Getting the DB name of the spatialite table loaded inside QGIS using PyQGIS

texts: I have a table loaded in QGIS from spatialite DB. When i see the layer properties it is dbname='/home/farook/Spatialite/Cities.sqlite' table="cities" sql=. Now i want the dbname alone (i.e) the entire path '/home/farook/Spatialite/Cities.sqlite' for whatever layer loaded inside QGIS dynamically using PyQGIS. Is there any options for the same?


texts: Shapefile buffer assigning to each buffer a different length

texts:
I have drawn a set of lines (corresponding to a network) in a shapefile. I would like to buffer them, but each line has a different buffer length, so I thought to do it with Python: for each branch, identified by a particular ID number, perform the buffer up to a certain distance that you can find in an attribute column.

To do so, I would like to use:
Buffer_analysis {all_my_branches, a_new_defined_layer, an_attribute_column_inserted_in_the_shapefile}

I imported arcpy, I set as arcpy.env.workspace my shapefile, but when I ask for the list of feature classes I get the answer:

"function ListFeatureClasses at 0x1948B770"

So I guess I have no features in my shapefile. How can I therefore do to perform the buffer I want? Can it be done on a shapefile?

I am new to ArcGIS and Python.


texts: How do I add projection information to this JPG and convert to geoTIFF in GDAL?

texts:
I have a .jpg image (weather satellite) that needs to have projection information applied to it and converted to geoTIFF in EPSG:3857, but due to the complexity of the projection I have no idea where to start using GDAL.
The .jpg (image.jpg) can be found here for reference, if it matters.
The projection of the image should be the following (taken from a netCDF file of the same data) geostationary projection:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Spheroid",6378137,298.2572221]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",-75],
    PARAMETER["satellite_height",35786023],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"]]

I have tried about 50 combinations of scripts using gdalwarp and even though I like to show effort instead of blatantly asking for help, doing so would only clutter my question. I am pretty lost.
Is this possible to do, or do I need more information than the gdalinfo I pasted above?
Edit: I also know what the bounding box should be : (after being reprojected to EPSG:3857)
West    -151.8860639
North    56.77489167
East    -52.91295
South    14.56403889

Edit 2: This is a netCDF file that has the same projection needed for the .jpg : https://s3.amazonaws.com/noaa-goes16/ABI-L1b-RadC/2020/211/05/OR_ABI-L1b-RadC-M6C13_G16_s20202110526162_e20202110528547_c20202110529074.nc


texts: Filter/sort image collection according to one of its band values in Google Earth Engine

texts:
I would like to filter/sort two sub image collections out of one image collection. The sub collections would be filtered/selected by higher 20% and lower 20% of one of its band values.

In other words, what I try to do is:

load image collection;
calculate and add NDWI index to this collection;
use that new ndwi index band to select/sort/filter two new image
collections.

One would consist of images with the lower 20% NDWI index values; and the other would consist of images with the higher 20% of NDWI index values.

So far I tried to use the arraySlice. It works, but unfortunately I only managed to calculate an image, which consists of the mean of the images with the lower/higher 20% ndwi index.

My next step would be to calculate, for either of these two new image collections, several other indices and apply a combined reducer. Therefore I need an image collection, and not just a composite.

Basically, I am just trying to use the ndwi index values as a filter for my image collection.

Please find a toy script here: https://code.earthengine.google.com/8101949c01e192cc8caba4dccff25583


texts: QGIS defaults new layers to not defined CRS (EPSG 7791)

texts:
Since I updated QGIS to 3.12 whenever I load a new layer into QGIS it sets the CRS to EPSG 7791. The same CRS is used as default when I want to export a layer, I always have to change the CRS from EPSG 7791 to my desired EPSG 3044 before I can export the layer. I never used EPSG 7791 before. I have set the Option under "Settings" -> "Options" -> "CRS"  to EPSG 3044. The project CRS is also 3044, so I really do not understand where  EPSG 7791 comes from. Anyone experienced the same issue?

Edit: Even after complete uninstall (and removing the remaining "profiles" folder) and reinstall the issue remains.


texts: How to narrow the result from Nominatim with only cities and countries, no additional areas with the same name?

texts:
I'm trying to get WKT or GeoJson polygons of cities and countries from Nominatim (for example, Paris and France). I send a GET request but it returns me several results named "Paris" or "France". How may I narrow results only to a capital and a country. Here's my request url:

http://nominatim.openstreetmap.org/search/fr/paris?format=json&polygon_geojson=1

And I get an array of 10 objects where only one is the capital:

class: "boundary", type: "administrative".

The rest are:

class: "place", type: "isolated_dwelling"
class: "place", type: "hamlet",
class: "place", type: "neighbourhood"
and so on...

How to select only the desired city and only the desired country ignoring rest of results with the same name?


texts: How to add a custom SymbolLayerV2 to a Vectorlayer in QGIS

texts:
I'm trying to add a symbolLayer with a line-decoration (line-arrow) to a vectorlayer of linestrings. what i am trying at the moment is to create a new layer and set it's QgsLineSymbol with a list of symbollayers.

# Symbology
        simple_props = { 'width' : '1', 'color' : '255,255,255' }
        deco_props = { 'width' : '0.26', 'color' : '0,0,255' }

        simple = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("SimpleLine").createSymbolLayer(simple_props)
        deco = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("LineDecoration").createSymbolLayer(deco_props)
        s = QgsLineSymbolV2([simple, deco])

        self.ll.setRendererV2( QgsSingleSymbolRendererV2( s ) )

Every time QGIS trys to set the QgsLineSymbolV2 it crashes with a win32-exception. what am i doing wrong, or is there a possibilty to just add a new symbollayer to a existing vectorlayer?


texts:
I have shapefile with few polygons. Some of them are self-intersecting. In ArcMap I can use "Repair geometry" tool, and it works well. But problem is, that for that task some users need to use open source or free tools.

I found tutorial on link here how to do that, but it is too hard to do it without explaining how it works.

Is there any other fast, one-click process for repairing geometry, like in ArcGIS software?


texts: How to change Leaflet marker color on mouseover?

texts:
I'm having a hard time figuring out why the Leaflet method setstyle will change the color of a polygon but not the color of my markers.

Polygon works fine:

But the markers won't change colors:

I want to be able to mouseover a marker and changes its color.  It seemed like setStyle would do this.  But I keep getting layer.setStyle is not a function

I'm using Angular and Leaflet together to make the map (I'm using the angular-leaflet-directive).

Here's the mouseover part of the code:

$scope.$on("leafletDirectiveMap.geojsonMouseover", function (ev, leafletEvent) {
        pointMouseover(leafletEvent);
    });

    function pointMouseover(leafletEvent) {
        var layer = leafletEvent.target;
        layer.setStyle({
            weight: 2,
            color: '#666',
            fillColor: 'white'
        });
    }

When the marker get's moused over, it fires pointMousever which then tries to invoke setStyle on the LeafletEvent.target.  I've console.logged the LeafletEvent and there is indeed the Target part of it:

Why would setStyle work for the polygon and not for the marker?  What's the way to change the marker's color?


texts: How print map in ol3 with top resolution?

texts:
I'm developing a GIS app using OpenLayers-3. In this app, users draw features and print the map.

Suppose the user draws his features in a large area. Now to print all his features, he would zoom out until all his features are in the screen. In this zoom level, quality of print will be very low.

I want the user to set a zoom level and then print all features in this zoom level. In this state, it's possible some features will not be visible. Finally the result should be a high quality image, with resolution greater than screen resolution.

How can I Implement that?


texts: Geoserver shape simplification performance

texts:
So, while wondering about shape simplification in the eternal struggle for optimization, I find out (here on SE) that Geoserver and PostGIS can and will (by default) perform the trick on the fly (see Small features don't draw when zoomed out using Geoserver WMS), all by themselves.

A number of questions come to mind (does it work with shapefiles too ? at which point of the process does it happen ? obviously the renderer must simplify the shape in any case, simply because of available pixels - but does a non-simplified shape that gets simplified on rendering have the same weight resource-wise as one that's pre-simplified ?), but to pick one:

how much of a resource hit is it, compared to preparing simplified versions of the geometries beforehand ?

The spider sense says it should be significant, but Geoserver has surprised me several times already in the brief time we have known each other.


texts: Publish Vector and Raster (GeoTiff) layer in Geoserver SLD

texts:
Base layer and overlay layer are same when opening in ArcMap.

When publishing the same layers in geoserver , I get a shift of 27 meter.
I have attached an example image to show the problem.


texts: ArcGIS drops some elevation values when converting to 3D or DWG format

texts:
I have a shapefile of 2' contours that I've converted into dwg format.  I have an existing elevation field (double) in the shapefile that's populated with my elevation values.  When it converts to dwg format there are some elevation values that don't carry over and populate as 0.  These records all have an Entity field of Group instead of 3DPolyline.  I've also tried the "Feature to 3D by Attribute" tool with the shapefile and it also assigns some of the contour lines to this Group entity and the elevation values don't carry over.  There are approximately 50 of my 650 contour lines that don't carry over correctly.  I've tried this with several shapefiles both in ArcGIS desktop 10.5 and ArcGIS Pro and the result is the same.

Why is this conversion dropping some of my elevation values and what can I do to correct this? Many of these dropped values are for a series of adjacent contour lines making it difficult to even extrapolate the missing values.


texts:
The Planet says that we need to add API key to the environment.
I am trying to do it as:

export PL_API_KEY=a3a64774d30c4749826b6be445489d3b

and getting SyntaxError: invalid syntax

What I am doing wrong?

I am typing it in cmd.exe.

I am trying to do:

pip install planet

It is not working either

Somehow, I feel that all these errors have the same source, but I can't fix it.


texts: How to make lines styles with offset join properly?

texts:
I want to draw the lanes of a two-lane street in different colors. I tried to apply an offset to a symbol layer which works fine for a single line but the ends of two connected lines do not join properly. Depending on the angle, the lines overlap or there is a gap between the lines.

Is there any way to solve my use case?


texts: Why new feature disappears when editing polygon?

texts:
I try to create internal polygon within existing polygon.
After I created a ring, I add a new feature inside the ring.
Suddenly, this new feature disappears. I confirmed that this feature went to the origin( 0,0).

The polygon and project file are used UTM as projection.

If somebody knows, why it is happening, let me know the reason and solution.

QGIS:1.8.0
OS:Windows 7


texts: Filtering out attributes fields of a vector layer using PyQGIS

texts:
I have a vector layer with 50 columns/fields. I have been trying to find a way of selecting 6 columns of particular interest (by name or index) and then creating a new layer from that selection, such that my new layer only has those 6 columns within it. I am hoping to achieve this with PyQGIS
At present, the only way I've solved this is to create a copy of the layer and then, following the documentation, delete attributes I don't want to keep from that copied layer. Seems very clunky. I would prefer a method that does this by selection rather than deletion

caps = layer.dataProvider().capabilities()
if caps & QgsVectorDataProvider.DeleteAttributes:
    res = layer.dataProvider().deleteAttributes([43,42])
    layer.updateFields()

Is there a more elegant way of doing this?


texts: I am trying to build a chloropleth map in Javascript with zoom functionality. Map will need to have a chloropleth layer of all countries. On zooming in it will need to show a chloropleth of the next adminstrative areas across all countries. Can you please give me pointers around how do to this. How long does loading all-states GeoJSON take and how do we go about doing this.


texts:
I'm new to working with GeoJSON and have been reviewing the ArcGIS page for AsShape, but the examples seem a little thin.

I'm currently pulling a WFS from a GeoServer instance in GeoJSON format.  I am using pycurl to do this as the service is https and it's easier to do authentication with pycurl.  I basically issue my url query and build a file (results.json).

I know how to easily build point feature classes from any text (including json) and build attributes using cursors, but the data I'm working with now has lines and polygons as well.

Could anyone give me some guidance, share a snippet, or something that basically shows me how to get into the .json file with AsShape, or uses AsShape successfully with the json.load command?


texts:
Does anyone know where I could download geoids? Specifically I need the Swiss Geoid CHGeo2004.

I found this link: https://opendata.swiss/de/dataset/geoidmodell-chgeo2004-in-ch19031

But it only includes a RGB TIFF and some contour lines, which I guess is only for visualization.


texts: How to maintain Shape_Length field afer Spatial Join of Point to Polyline feature classes?

texts:
I need to get for each point the length of the line that joining to the point.

How to maintain the "Shape_Length" field of polylines when I'm Spatially Joining point to lines feature classes ? When I'm using the normal "Spatial join" tool this field isn't appear at the point output layer.
 Apparently because this is the default field of line feature class it disappeared at point layer.


texts: I currently experiment with an affordable L1 GPS RTK system based on u-blox NEO-M8T with RTKLIB and higher quality patch antenna (with ground plane). I wonder what is the role of L2 signal in an RTK system (since the expensive RTK setups seem to use L1+L2) - how much precision/accuracy will it add to the system?


texts:
I want to add new column to attribute table of vector layer via Python console. I'm trying to do this following instructions from http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html.
But while defying caps

caps = layer.dataProvider().capabilities()

I'm getting error:

  AttributeError: 'NoneType' object has no attribute 'dataProvider'

I've found some solutions referring to setting up QGIS Python environment correctly. (Getting dataProvider from vector layer outside QGIS). But this solution crash my QGIS.

Any idea what should I do to access (?) (or rather 'get?') dataProvider?


texts: Intersecting overlapping polygons individually in ArcGIS Desktop

texts:
Map Layers:

Census Tract ShapeFile with Population Data
Set of Store Locations with 5 Mile Buffers (Significant Overlap Between Buffers)

I have create feature layers for both and I am trying to intersect the population data and the buffer data.  However, the overlap between the buffer polygons is leading to significant duplication/over counting, as the polygons draw from the census tract layer as well as the overlap with other polygons.

Is there a way to intersect the polygons individually so that they are only drawing from the census tract data?


texts:
I have a SQLlite database containing cached tiles with X,Y,Z fields and I want to use it in MapCache.
My Y field  is appearing upside down in Google tile index but I want to inverse it like this:

http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/

I'm searching for a way to use my databse in Mapcache or inverse Y field in database to adapt tile indexes with mapcache tile index.


texts:
I have GridCoverage2D object built from the two-dimensional float array:

float[][] myData = new float[dataH][dataW];
...
// fill myDdata array
...

ReferencedEnvelope env = new ReferencedEnvelope(-180, 180, -90, 90, DefaultGeographicCRS.WGS84);
GridCoverage2D gc2d = gcf.create("coverage", myData, env);

I have a number of no-data points in myData array. No data value is -999.
Is it possible to interpolate myData and "fill" these no-data values using Bilinear or Bicubic interpolation?

I have tried:

javax.media.jai.Interpolation interp = Interpolation.getInstance(Interpolation.INTERP_BILINEAR);
GridCoverage2D gc2d_interpolated = (GridCoverage2D) Operations.DEFAULT.interpolate(gc2d, interp);

However with no results. What I am doing wrong? Should I set NO_DATA some way to my interpolation?

******* EDITED *******

It is difficult to reproduce full section of code that's why I prepared some example reproducing the problem.
I attached small subset of my data as a GeoTiff image

It is my example:

File inputFile = new File("E:\\tmp\\img\\subset1.tiff");
File output1File = new File("E:\\tmp\\img\\no_interpolation.tiff");
File output2File = new File("E:\\tmp\\img\\bilinear.tiff");

AbstractGridFormat format = new GeoTiffFormat();
GeoTiffReader reader = new GeoTiffReader(inputFile);
GridCoverage2D c = reader.read(null);

GeoTiffWriter w = new GeoTiffWriter(output1File);
GeoTiffWriter wi = new GeoTiffWriter(output2File);
w.write(c, null);

javax.media.jai.Interpolation interpolation = Interpolation.getInstance(Interpolation.INTERP_BILINEAR);
GridCoverage2D ci = (GridCoverage2D) Operations.DEFAULT.interpolate(c, interpolation);

wi.write(ci,null);

w.dispose();
wi.dispose();

Both resulting images with and without interpolation are attached. There are no difference between them.


texts:
My vector data is based on EGM2008 and I want to add local geoid to arcgis and project vector data from EGM2008 to local geoid.

If it is not possible to add custom local geoid in arcgis, how can I project/ transform 3d vector data to local geoid?


texts: Ship detection with ArcGIS Pro - deep learning

texts:
Currently I try to develop and apply a deep learning model to detect vessels. I imported aerial imagery (3 bands, Color Infrared, 25 cm resolution, ECW files) and created manually about 200 training samples of ships. As software I used ArcGIS Pro 2.5 in combination with Jupyter Notebook.

I exported the training samples to image chips and trained a model with 78% accuracy in Jupyter Notebook which looks promising.
As a last step I used the 'Detect Objects for Deep Learning' tool in ArcGIS Pro to apply my model. Part of the results are shown in the image below:

Did I something wrong that I get different predictions for actually one individual ship?
I hoped to get something like this:

Any suggestions how could I improve my results?


texts: ArcGIS Runtime .NET does not make a GetLegendGraphic request

texts:
I'm using ArcGIS runtime SDK .NET to populate maps with WMS Layers

This is a bit strange. If I load a WMSLayer on the map without specifying any layer, I see how it starts to recover all the legends of the server, one by one. On the contrary, if I specify which layers I need (setting field Layer) I see that it adds the layer to the map, but does not make the request to read the legend

If Layers object is empty:

Do a GetCapabilities request
Do one GetLegendGraphic request per layer
Do a GetMap request

But if I set Layers object

Do a GetCapabilities request
Do a GetMap request

In both cases ShowLegend is true


texts: Where can I find paved/unpaved road surface data for the US?

texts:
I'm currently working on a small program that evaluates OpenStreetMap data to calculate how twisty or curved roads are in order to generate maps of roads that may be of interest to motorcycle riders or other driving enthusiasts: https://github.com/adamfranco/curvature/wiki

The major limitation I've found is a dearth of road surface data that indicates whether roads are asphalt, concrete, gravel, or at a minimum paved/unpaved. Are there any non-proprietary data sets that contain this information for either the whole US or particular states? The TIGER dataset doesn't seem to contain road-surface information as far as I can tell.


texts:
I have loaded some rasters into Postgis. I did a lot of ST_XXXX (scale, snap to grid,resample) magic to make them aligned. After executing:

select
r1.rid,r2.rid
,ST_NotSameAlignmentReason(r1.rast,r2.rast)
,ST_SameAlignment(r1.rast,r2.rast)
from
  rasters.r r1,
  rasters.r r2
where not ST_SameAlignment(r1.rast,r2.rast)

I got nothing, if I remove:

  where not ST_SameAlignment(r1.rast,r2.rast)

all I got is

  "The rasters are aligned"

Which seems ok, but when I execute:

SELECT ST_Union(rast,'FIRST') from rasters.r

I got:

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

Why ?

Postgres 9.4

"POSTGIS="2.1.8 r13780" GEOS="3.5.0-CAPI-1.9.0 r4090" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER"


texts: Publishing mosiac datasets for ArcReader?

texts:
Has anyone had luck with publishing their own imagery into ArcReader 10.2 in a mosaic dataset?

I have had success publishing the my own imagery in tiles but not in a mosaic form.

Does anyone have any suggestions or tips?


texts: QGIS: can I somehow display a virtual raster of 7000 zipped GeoTiffs?

texts:
I've built a virtual raster of 7000 zipped geotiffs (each about 3600x2400 px). The GeoTIFFs have overviews.

I've tried to open it with QGIS desktop application (via Add Raster Layer).

Alas, QGIS hung up.

Can I somehow show such a layer? Maybe I should split it to tiles?


texts: Looping spatial join in arcpy - order problems

texts:
I developed a python script where the user inputs

a polygon shapefile
multiple point shapefiles

and within the script the polygon shapefile is divided into 1000mx1000m grid cells and the point shapefiles are joined to the polygon where the intersecting points per cell are counted.

The spatial join is done within a for loop and after the joins are finished, all fields that are not required and do not start with "Join" (like Join_Count, Join_Cou_1, ...) are deleted. Then, the output shapefile is stored.

I tested the code and everything was fine, but now I detected a huge issue. The order of the Join_Count, Join_Cou_1, Join_Cou_2, etc. does not correspond to the order in which the point shapefile were entered.

Below the input order and output order can be seen. I tried it multiple times, but it seems there is no logic behind the order, so not alphabetical, Input order, reverse order, or anything. It seems like it would be ordered randomly.

Bank -> Join_Count (Bank)

Kino -> Join_Cou_1 (Post)

Park -> Join_Cou_2 (Park)

Post -> Join_Cou_3 (Kino)

Why is this happening? Shouldn't it be in the same order as the input? Of course I cannot apply my code then, because the user would have to look manually which of the inputs are at which position.

Here is my code:

import arcpy

studyarea = arcpy.GetParameterAsText(0)
inputPoints = arcpy.GetParameterAsText(1)
outputLocation = arcpy.GetParameterAsText(2)

arcpy.GridIndexFeatures_cartography('in_memory\studyarea_grid', studyarea, 'INTERSECTFEATURE', 'NO_USEPAGEUNIT', "#", 1000, 1000)
studyAreaGrid = arcpy.Clip_analysis('in_memory\studyarea_grid', studyarea, arcpy.Polygon)

inputPoints_set = inputPoints.split(';')

spatialJoin = studyAreaGrid
for pointShp in inputPoints_set:
    spatialJoin = arcpy.SpatialJoin_analysis(spatialJoin, pointShp, arcpy.Polygon, 'JOIN_ONE_TO_ONE', 'KEEP_ALL', '', 'INTERSECT')

fieldNameList = []
fieldObjList = arcpy.ListFields(spatialJoin)
for field in fieldObjList:
    if field.required == False:
        if not field.name.startswith('Join'):
            fieldNameList.append(field.name)

arcpy.DeleteField_management(spatialJoin, fieldNameList)

arcpy.CopyFeatures_management(spatialJoin, outputLocation)

I again tried to use fieldmappings, see code below.
It works, meaning that only the Join_Count fields are passed on, but also the order is not the same as the Input order.
Why does this happen? Should not fieldmappings work??

for pointShp in inputPoints_set:
    fieldmappings = arcpy.FieldMappings()
    fieldmappings.addTable(spatialJoin)
    fieldmappings.addTable(pointShp)
    for field in fieldmappings.fields:
        if field.required == False:
            if not field.name.startswith('Join'):
                fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex(field.name))
    spatialJoin = arcpy.SpatialJoin_analysis(spatialJoin, pointShp, arcpy.Polygon, 'JOIN_ONE_TO_ONE', 'KEEP_ALL', fieldmappings, 'INTERSECT')


texts: Import database (Excel file) into PostgreSQL

texts:
In order to import my database (excel file) into PostgreSQL, I converted the Excel file into ".csv" file. Then, I used the application: "Import Data For PostgreSQL" and  I could import my database but only the first 100 lines; knowing my database contains 80753 registration.

I tried  to import the database by SQL query, but an error message appear and I don't know how to solve it:

  ERREUR:  séquence d'octets invalide pour l'encodage « UTF8 » : 0xcce5.


texts: Converting text to proper case using ArcGIS field calculator?

texts:
I am working on a set of data that has over 40,000 records on ArcGIS and I want to use a Proper Case for several fields/columns on the table.

Is there a Python Script that can help me automate this process as I desire to save time while working on this records?


texts:
I have generated the following command line in ERDAS, and want to run it in MS-DOS.
Getting an error- "system cannot find the path specified".
I have checked the path by opening it in windows explorer, it navigates to the corresponding file,
have tried multiple options, by changing the location or going through .bat approach, still it won't work!

C:/ProgramFiles/Intergraph/ERDAS IMAGINE 2014/bin/x64URelease/modeler.exe
@"C:/Users/WRD/AppData/Local/Temp/ImagineCommandLine_modeler_3780003780"


texts: QGIS server 2.6.1 upgrade to 2.8.1 error

texts:
I plan on upgrading QGIS server from 2.6.1 to 2.8.1. At the time of 2.6.1 is used, however, to 2.8.1 found the map service can not be used, the emergence of internal server error 500.
C:\OSGeo4W64\apps\qgis\bin\qgis_mapserv.fcgi.exe - FastCGI Process unexpectedly quits

I will qgis server as a fcgi module placed on iis, rather than on Apache. When I use QGIS 2.6.1 version, my service request is normal. this is the way my request:

http: //localhost/qgis_mapserv.slg SERVICE = WFS & version = 1.0.0 & request = getfeature & map = C:? /gis/zhongguo.qgs&typename=zhongguo
This is my request to return an error

Modules
FastCgiModule

Notification
ExecuteRequestHandler

Handler
qgisserver

Error Codes
0x00000001

Seems to be the application initialization error

Who can help me?


texts: Getting error when I add drawing feature to the map

texts:
I use OpenLayers 6.
I try to add to the map interaction a class to draw features on the map.

Here is MapDraw class:

function MapDraw() {

        this.vector = new ol.layer.Vector({
            source: new ol.source.Vector()
        });

        Draw.call(this, {
            source: this.vector.getSource(),
            type: "Point"
        });

        this.on = Draw.prototype.on;
        this.on('drawend', function (e) {
            eventsHandlers.featureComplete;
        });
    }

    if (Draw) MapDraw.__proto__ = Draw;
    MapDraw.prototype = Object.create(Draw && Draw.prototype);
    MapDraw.prototype.constructor = Draw;

    $("#popupShapes").popup("close");

    return MapDraw;
}(ol.interaction.Draw));

Here how I add class above to the map interaction:

map.addInteraction(new MapDraw());

But when add interaction occurred on this block:

    Draw.call(this, {
            source: this.vector.getSource(),
            type: "Point"
        });

I get the following error:

Uncaught TypeError: r.setActive is not a function

Any idea why i get the error and how to fix it?


texts: How to update coordinates by values in columns?

texts:
I try to make a script which changes the coordinates given in format DMS (like N54 24 04.50, E018 34 17.20 to 18.6889, 54.3320) to DD format which is suitable for working in QGIS. The conversion is pretty easy but what I want to really do is to change a point position basing on coordinates given in columns. Like here:

So I want to create a feature, then type in LAT and LON columns the desired position, columns 'x' and 'y' shall take these coordinates and convert them to DD format. And the last step is to change the feature position basing on columns 'x' and 'y'.

I know that I can retype the coordinates or recalculate them in csv. sheet and load again to QGIS but it's time consuming. So is there a possibility to 'set' a columns from which the feature will take the $x and $y of geometry directly in QGIS?


texts: Polygon to raster with extent equal to the raster one?

texts: I need to produce a range of raster layers with the same extents to be able to apply "Weighted overlay". The problem is when converting from polygon to raster the output raster has 1 value for the polygon but Nodata value for the rest of the extent. How could I manage to get a data value instead for this area?


texts: Creating new layer from selection in QGIS?

texts:
Is it possible to create layers from selection in QGIS?

I am using QGIS 1.7.4-Wroclaw.

To make it easier for others to see this, here are some pics. I used QGIS 2.0.1-Dufour.

SELECT

SAVE

ADD


texts:
I have a feature dataset with 317 points (UpperPoints) and another with 394 (LowerPoints) points. Now I want to calculate the cost paths as in this example here. However there are multiple files used and i just want to iterate trough my shapefiles. In the end I want the cost paths from every single UpperPoint to every single LowerPoint. I know that I have to calculate the cosDistance rasters for every lower Point which will need a lot of computation time, but that's not an issue. My main problem is how to extract the points from the shapefiles and iterate trough these to create the CostDistance rasters. I assume that my solution, just using the for loop is too simplistic...

PointsHIGH = "C:\\Dogon\\lcp.gdb\\Pts_rangeH_high"
PointsLOW = "C:\\Dogon\\lcp.gdb\\Pts_rangeH_low"
costsurface = ("C:\\Dogon\\lcp.gdb\\Slope_asc")

outputCostDist = "C:\\Dogon\\lcp.gdb\\CstDstLOW"
outputCostBackLink = "C:\\Dogon\\lcp.gdb\\BckLnkLOW"
outputCostPath = "C:\\Dogon\\lcp.gdb\\CostPath"

for points in PointsLOW
    outputCostDist = arcpy.gp.CostDistance_sa(points, costsurface, outputCostBackLink)
    outputCostPath = arcpy.gp.CostPath(points, outputCostDist, outputCostBackLink,"BEST_SINGLE")

Further I want an output of all CostDistance rasters and CostPaths, how could I do this?


texts:
I am using QGIS 2.18 and I am wondering how to find the NW corner of a polygon and create an attribute of that point number. See the attached image.
For example, please node 8 is the NW corner for Grid #1, 13 is the NW corner for Grid 5 and so on. any help with this would be wonderful. The two files I am working with is

Grid = 'grid'
and

Vertex = 'nodes' (Bold)


texts: Setting "Disable tiled raster layer exports" option to export to PDF using PyQGIS

texts:
I need to use the setting Disable tiled raster layer exports from the GUI in my Python code, but I can't seem to find it.

exporter = QgsLayoutExporter(layout)
exporter.exportToPdf(path, QgsLayoutExporter.PdfExportSettings())

How can I do this?


texts:
I'm still a js noob.

I am currently developing a map using the Mapbox GL JS API with their geocoding plugin, however, I'm am not seeing any examples of adding a marker after the return.

I did find a previous thread (Mapbox Geocoding and Markers) with the same question, however, this is referring to the js api and appears to be using a different syntax.

Here is my code thus far:

<!DOCTYPE html>

<meta charset='utf-8' />
<title></title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.12.2/mapbox-gl.css' rel='stylesheet' />

<link rel='stylesheet' href='https://api.mapbox.com/mapbox-gl-js/plugins/mapbox-gl-geocoder/v0.0.0/mapbox-gl-geocoder.css' type='text/css' />

Here is the link to the API: https://www.mapbox.com/mapbox-gl-js/api/


texts: How does OGR detect what GEOS version I have?

texts:
I'm attempting to solve this error:

ERROR 6: GEOS 3.8 or later needed for MakeValid.

The version of GDAL I am running is 3.0.4 which was downloaded as a .whl from https://www.lfd.uci.edu/~gohlke/pythonlibs/

The version of GEOS that ships with this .whl is unknown so I assume it's <3.8.

To that end I have downloaded GEOS 3.8.0 separately and managed to compile it sucessfully and generate the geos_c.dll and geos.dll to replace the ones in my current python installation.

The only issue is that the error remains still remains when trying to run the code.


texts: Creating map book series with blank pages between maps using ArcGIS Pro?

texts:
I am looking to create a series of maps each one highlighting the camping pitch,. I would like to put the number of the pitch highlighted on either page before or after to allow duplex printing of map on one side and pitch number on the other side?.

Is this possible using ArcGIS Pro with or without python?

See image of camping pitches -


texts: Bulk loading Picture Marker Symbols into .style for ArcMap?

texts:
I have close to 1000 PNG files that I want to make Picture Marker Symbols stored in a .style file.  I am very familiar with Microsoft Access, Python and ArcGIS and I know some C# and .NET.  When I look at the symbols in the .style with Access, they are stored as a Long binary data.

Does anyone know how I can get a PNG into the .style file without loading each one individually?


texts: Change date format into readable format in GEE

texts:
I have code in GEE which create table with the sentinel2 mean NDVI data for certain period.
my problem is that the table in the end has two column with the dates, one is the name of each image (not easy to read- 20191004T051651_20191004T052839_T43QGC_00000000000000000000 for example) and the other has also dates but not readable ( for example : 1573623201083)

I have tried to change the format in different places in my code so it will get the correct time values, e.g:

// Collect region, image, value triplets.
var triplets = filter.map(function(image) {
  return image.select('NDVI').reduceRegions({
    collection: geometry.select(['Name']),
    reducer: ee.Reducer.mean(),
    scale: 20
  }).filter(ee.Filter.neq('mean', null))
    .map(function(f) {
      return f.set('time', image.get(ee.Date('system:time_start')));
    });
}).flatten();

print('triplets', triplets);

  FeatureCollection (Error) Date: Bad date/time 'system:time_start'

and:

      // Get the list of all features with a unique row ID.
      var values = ee.List(row.get('matches'))
        // Map a function over the list of rows to return a list of
        // column ID and value.
        .map(function(feature) {
          feature = ee.Feature(feature);
          return [feature.get(Name), feature.get('mean')];
        });
      // Return the row with its ID property and properties for
      // all matching columns IDs storing the output of the reducer.
      // The Dictionary constructor is using a list of key, value pairs.
      return row.select(ee.Date([time])).set(ee.Dictionary(values.flatten()));
    });
};

  Invalid argument specified for ee.Date(): time

This is link to the script:
https://code.earthengine.google.com/f4bb8bcac7d2e351523b33d9401a8618

it has already the geometry , in case it doesn't work let me know

My end goal: to have readable date column in my table


texts: Get the geometry from the text in a label

texts:
I need to be able to convert the outline of label text to a polygon shapefile so I can use it with ArcGIS Online.

Can I somehow get the geometry of the text values using pyqgis?


texts: Show google coordinates in leafLet with WMS Geoserver

texts:
I am new to GIS and have a problem with using geo-coordinates which I got from  Google. I have some lat, lng data in my database(MSSQL) which I obtained from Google API (I think it should be in 3857 format), however, I needed to create a layer in GeoServer, so I created a view as

   geometry::STGeomFromText('Point(' + CAST(dp.lat AS nvarchar(50)) + ' ' + CAST(dp.lng AS nvarchar(50)) + ')', 3857) AS Expr1

In Geoserver, when adding a layer I use the following setting which I believe there is something wrong with it and I am not sure which Compute shall I click and in what order.

In the publishing, I also select point as default style.

and finally in leafLet

 var mywms = L.tileLayer.wms("http://localhost:9090/geoserver/localhost/wms", {
        layers: 'localhost:DP_View',
        format: 'image/png',
        transparent: true,
        version: '1.1.0'
         ,SrsName: 'EPSG:3857'
    });
    mywms.addTo(map);

They either sit in the middle or on top of the map.


texts: Error Reporting Registry Location of ArcGIS for Desktop?

texts:
I need to stop ArcMap from error reporting during an Arcpy loop. The error report dialog stops my loop and I need for it to go away, permanently

I found this link which references the registry location as 'HKCU\Software\ESRI\Settings\ErrorReports'.

The actual folder hierarchy in my regedit is much different and I cannot locate the ErrorReports registry

What is the actual location of the registry responsible for error reporting in ArcMap 10.2.2?


texts: How to find point within polygon with in TURF


texts: Is there a "Select by Attributes" plugin/widget for ArcGIS Server JS API?

texts:
I'm hoping to find some existing code which emulates the "Select by Attributes" dialog in ArcMap, within the ArcGIS Server JavaScript API.

Writing a hard-coded QueryTask is easy - I'm looking for some code which will automatically handle any layer on the map, forming a Query on-the-fly using any attributes, ideally returning the results into a datagrid.

Clever people have written Table of Contents and Identify widgets which handle any layers in the map. Has anyone written something similar for "select by attributes"?


texts: When to use ModelBuilder over Python scripting and vice versa?

texts:
I am new to Python scripting in ArcGIS but not new to ModelBuilder.

I would like to know the benefits of Python vs ModelBuilder.

When should we write Python script for GIS automation rather than use ModelBuilder?

What are the automation capabilities of Python that we cannot find in ModelBuilder?


texts: How to start QGIS after building & install on ubuntu 12.04 64bit?

texts:
I am a newbie here.

I am running ubuntu 12.04 64bit, and build after the instruction, but i don't know  why can not launch qgis. Seem i successfully build $ install it.

below is the linux shell:

jack@jack:~/apps/bin$ sudo su
root@jack:/home/jack/apps/bin# qgis
The program 'qgis' can be found in the following packages:
 * qgis
 * qgis-plugin-grass
Try: apt-get install <selected package>
root@jack:/home/jack/apps/bin# ls
qbrowser  qgis  qgis_bench
root@jack:/home/jack/apps/bin#


texts: Transfering color style from one raster to other rasters using QGIS?

texts:
I have a raster set and I need to classify them using the "pseudo color singleband" with "RdYlGn". What I'm trying to do is classify a raster and copy the style to the others, but it does not work because it limits the range to the values of the first raster I used as a template.

Is there a way to classify the raster, copy the style to others and it refresh the range automatically?


texts: Ask ArcMap to recognize updated CSV file content

texts:
I have a CSV file that I am updating outside of ArcMap (v10.4.1).
How can I get ArcMap to see the revised content instead of the old content?
ArcMap does not automatically see the updated version.
My current solution is to exit ArcMap and start a new session, which is time consuming.
I have already tried:

removing the CSV from the Table of Contents and reintroducing using drag-and-drop from the Arc Catalog pane,

removing from TOC, then reintroducing using the "Add Data" tool ("Standard" toolbar).

starting a new MXD using the "New" tool ("Standard" toolbar).

converting to DBF by right-clicking on the CSV in the Arc Catalog pane.

This problem seems  to be non-reproducible. I am not able to reproduce minutes after the original post.
I am running on a network which has noticible interruptions a few times a day, typically. I wonder if ArcMap will uses a cached copy if it cannot reload from the network.


texts: Create Python script that reads points(X,Y) from a fc, moves them in a user-defined distance and saves them in a NEW fc

texts:
i like to create a standalone Python script that reads points(X,Y) from a fc, moves them in a user-defined distance and saves them in a new fc. I tried to do the steps with the data access module, but i dont know so far how to go on, especially with the part "moving the XY-points" .

import arcpy
arcpy.env.overwriteOutput= 1
arcpy.env.workspace = "workspace"
outputFC = "New"
inFC = "inFC"

arcpy.CreateFeatureclass_management("workspace","New","POINT")
cur = arcpy.da.InsertCursor("New", "SHAPE@XY")

point=arcpy.Point()
cur2 =arcpy.da.SearchCursor("inFC",["SHAPE@XY"])
for row in cur2:
     print row
outputFC=row.GetValue("Shape")
point=outputFC.getPart()
point.Y = pnt.Y
point.X = pnt.X
array.add(point)
feat.shape=array
del row
del cur
del ins


texts:
I have little experience in scripting and I'm looking for a more efficient way  to replace the source file path of over 1000 shape files.

Is there a way I can do this without manually changing each file path?

As an example: I need to replace the S1 to S2 in the following file path

H:\GQRA screening tool\Shapefiles\S1 SOLID Ali _C10-12 0.5 m

this is the script i'm using, however the bit im struggling with is how to replace a partial string section of the file path.

import arcpy
... mxd = arcpy.mapping.MapDocument(r"H:\GQRA screening
tool\Arcmap\SwecoTemplate_SI_A3_S2.mxd")
... for lyr in arcpy.mapping.ListLayers(mxd):
...     if lyr.supports("DATASOURCE"):
...         if lyr.dataSource == r"H:\GQRA screening tool\Shapefiles":
...             lyr.findAndReplaceWorkspacePath(r"S1", r"S2")
... mxd.saveACopy(r"H:\GQRA screening tool\Arcmap\SwecoTemplate_SI_A3_New")
... del mxd

This is the updated script, however still not got it right yet.

import arcpy
... mxd = arcpy.mapping.MapDocument(r"H:\GQRA screening
tool\Arcmap\SwecoTemplate_SI_A3_S2.mxd")
... for lyr in arcpy.mapping.ListLayers(mxd):
...     if lyr.supports("DATASOURCE"):
...         if r"H:\GQRA screening tool\Shapefiles":
...             lyr.findAndReplaceWorkspacePath(r"S1", r"S2")
... mxd.saveACopy(r"H:\GQRA screening tool\Arcmap\SwecoTemplate_SI_A3_New")
... del mxd


texts: How to detect cloud's shadows with Sentinel2/SentinelHub?

texts:
I'm trying to detect clouds and their shadows on Sentinel2 images with SentinelHub service.

For clouds detection, SentinelHub team has developed a good tool but I would like to detect cloud's shadow too.

For the moment I use B08 band (NIR) with a threshold of 0.15, it's an empirical method I created. But maybe it exists a better way... Somebody have an idea?


texts: Merging features with the same values within feature class

texts:
I have a feature class of land parcels that contain "ownername" and "land_use_code" fields. Some of these features have the same values despite being two different parcels in different locations. For example, two parcels have ownername 'Jon Doe' and land_use_code '1001' but are on opposite sides of a town.
I want to merge those two parcels into a single feature based on those values and any others like that while also keeping all the other fields that are not relevant to the merge. Is there a way to do this for thousands of features at once in ArcGIS Pro?


texts: OpenLayers vector layer from GeoServer does not appear completely

texts:
I have successfully loaded the vector layer but the webpage only shows parts of it, some are not displayed. Any idea what might be the cause?

    var jsonSource_testLayer_1 = new ol.source.Vector({
    url: 'http://localhost:8080/geoserver/Testing1/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Testing1%3AtestLayer&maxFeatures=50&outputFormat=application%2Fjson',
    format: new ol.format.GeoJSON(),
});

var lyr_testLayer_1 = new ol.layer.Vector({
    source: jsonSource_testLayer_1,

    style: style_testLayer_1,
    title: '<img src="styles/legend/testLayer_1.png" /> testLayer'

});


texts: Some features listed in attribute table not visible on map in QGIS

texts:
My vector layer is a track file with points every 2-3 seconds. When added they all appear in the attribute table but not on the map view. There doesn't seem to be a pattern or any attribute that distinguishes those that are shown from those that aren't.

I'm using QGIS version 3.4.3-Madeira.

Any suggestions where I could find a solution?


texts: Earth Engine cast a string to a number for date filtering

texts:
I have a feature collection with the year of the campaign as a string "09/10".
To search through the require images I need to filter by dates turning the campaign into a year.

However I keep getting a message that complains I'm trying to do number operations on a string.

Number.int, argument 'input': Invalid type. Expected: Number. Actual: String.

The following script is heavily inspired by one of the answers given by Rodrigo E. Principe.

var mapfunc = function(feat) {

   // get the geometry of the individual feature
   var geom = feat.geometry();

   var myCampaign = ee.String(feat.get(campaign));

   var myYear = myCampaign.slice(0,2);
   // take the campaign of the feature and iterate through the correct
   // set of images
   var yearToIterate = ee.Number(myYear).toInt().add(2000);
   var subsetQKM = qKM
   .filterDate(ee.Date.fromYMD(yearToIterate,10,1),
   ee.Date.fromYMD(yearToIterate.add(1),5,25));

   // function to extract the value of the pixel for a feature
   var addProp = function(img, f) {

      var newf = ee.Feature(f);
      var date = img.date().format("Y-D");
      var b1 = ee.Number(img.reduceRegion(ee.Reducer.first(), geom, 250).get('sur_refl_b01'));
      var b2 = ee.Number(img.reduceRegion(ee.Reducer.first(), geom, 250).get('sur_refl_b02'));
      // Calculate the WDRVI the hard way
      var wdrviCoeff = ee.Number(0.2);
      var wDRVI = b2.multiply(wdrviCoeff).subtract(b1)
           .divide(b2.multiply(wdrviCoeff).add(b1));

      var wdrviDate = ee.String("wdrvi_").cat(date);
      //var b2Date = ee.String("sur_refl_b02_").cat(date);
      return ee.Feature(newf.set(wdrviDate, ee.Number(wDRVI)));
   }
   var newfeat = ee.Feature(subsetQKM.iterate(addProp, feat));
   return newfeat;
};

// add the data about the
var newft = ft.map(mapfunc);


texts: point in polygon starting with csv and shp in python

texts:
I have a csv file with two sets of lat/lon points - the first set of points is represented by columns V1(lon) and V2(lat),  and the second set of points by V3(lon) and V4(lat) (the csv does not have headers).

> head(merge)
                 V1               V2                V3               V4   V5     V6
1 -122.030277484219 38.0534546969437 -122.263688968328 38.0180951886275    0 Sunday
2 -122.012178329029 38.5883736394547 -121.269556181423 37.9550435479333    1 Sunday
3 -121.817242008376 38.8419619506953 -121.305704278949 37.2665873119052   10 Sunday
4 -122.100160967342 38.7339037710083 -121.839444904779 37.8529730771058  100 Sunday
5 -122.424897962683 37.3465530990606 -122.193097527866 37.2373493348652 1000 Sunday
6 -122.262026449704 37.6901105931389 -121.913111829827 37.0947943964057  101 Sunday

In python, I need to check the location of each point - and create a separate column in my .csv file containing the name of the feature in which the point is located. I have a shapefile containing polygons - lets call it counties.shp - so basically I need to find which polygon each point falls within.

I understand that first I have to import my shapefile, transform it to WGS84, import the csv, create points from lat lons, then loop over the points to search for overlapping features, although I am not sure about the syntax in python.  I have found the following resource but it needs to be modified.  Any ideas where to start?

#!/usr/bin/python
import ogr
from IPython import embed
import sys

drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file
ds_in = drv.Open("counties.shp")    #Get the contents of the shape file
lyr_in = ds_in.GetLayer(0)    #Get the shape file's first layer

#Put the title of the field you are interested in here
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm")

#If the latitude/longitude we're going to use is not in the projection
#of the shapefile, then we will get erroneous results.
#The following assumes that the latitude longitude is in WGS84
#This is identified by the number "4326", as in "EPSG:4326"
#We will create a transformation between this and the shapefile's
#project, whatever it may be
geo_ref = lyr_in.GetSpatialRef()
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)
ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)

def check(lon, lat):
    #Transform incoming longitude/latitude to the shapefile's projection
    [lon,lat,z]=ctran.TransformPoint(lon,lat)

    #Create a point
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    #Set up a spatial filter such that the only features we see when we
    #loop through "lyr_in" are those which overlap the point defined above
    lyr_in.SetSpatialFilter(pt)

    #Loop through the overlapped features and display the field of interest
    for feat_in in lyr_in:
        print lon, lat, feat_in.GetFieldAsString(idx_reg)

#Take command-line input and do all this
check(float(sys.argv[1]),float(sys.argv[2]))
#check(-95,47)


texts: How to Convert Number (Type Object) Into Type Number in Google Earth Engine?

texts:
I am creating a time-series chart of cloud cover. I calculate the cloud score for image using a custom function and store it as a cloud property.
Finally, I try to create a chart as follows:

Get the list comprised of each image's cloud score using aggregate_array

var cs_list = image_col.aggregate_array('cloud);

Finally, creating a series chart

var chrt = ui.Chart.array.values(cs_list, 1);
print(chrt);

And following is the script code:
var img_col = ee.ImageCollection("COPERNICUS/S2")
                    .filterDate('2018-05-01', '2018-07-30')
                    .filterBounds(region_of_interest_geometry);

var col_with_cloudiness = img_col.map(function(image) {
      var cloud_score = compute_cloud_score(image); // function to get cloud_score e.g 0.3454

      return image.set('cloud', cloud_score); // set new property
});

var cs_list = col_with_cloudiness.aggregate_array('cloud');
print(cs_list);  // JSON [0.2564 (type Obj), 0.0573 (type Obj), ...]. Want Array [0.2564 (type Number), 0.0573 (type Number)]

print(ui.Chart.array.values(cs_list, 1));

On running the code, I ran into the following issue: the cloud-score list contains numbers which are of type object. And so instead of getting line plot, several lines' legends representing each list element(score) gets created
How to convert them to numbers or create an array of cloud scores as actual numbers?


texts: Can a QGIS project (.qgs) be uploaded to work in map viewer in GIS Cloud?

texts:
Title says it all. Can a QGIS project be uploaded onto GIS CLoud to display all the properties i have asked it to display in QGIS?

if not is there an equvilant GIS Server i can upload previous QGIS project onto so they look like the look on QGIS?


texts: Adjusting heights within DEM raster for visualisation?

texts:
Problem
I have a bathymetry DEM which I want to display in 3D (ArcScene or Fledermaus).  The data set is one single raster but looks like it is a mosaic of several smaller datasets.  The problem is the absolute height values between these sub-sets varies and is quite prominent: resulting in the cliff like features in the below picture.

Goal
I want to 'correct' each of the subsets to create a more smooth and realistic scene.  I'm not bothered by the accuracy of the depths, this is just a relative correction for visualisation.
My Ideas
Originally I had hoped we had the individual sub-sets as separate raster and I was going to start from the largest sub-set raster, take an adjoining raster and add or takeaway the mean difference in the heights along their adjoining edge.  Maybe looking a bit at the difference in error from one side of the edge to the other.  For an indication the bin size is 100m and depths only change by 1-3m between cells within a sub-set.  Between sub-sets the difference is 30-100m.
Keeping his idea in mind, I then tried to separate the sub-sets back out ('un-mosaic').  Manually digitising will work but it's a lot of manual labour - too much to justify doing.  Picking out the edge has to be quite precise because if I am a cell out on either side the cliff edges will still be there.
I have been trying to find more automated ways of doing this and tried a lot of edge finding filters (3x3 and 5x5 using various values) and sharpening techniques although my experience of working with rasters in this way is quite limited.  It sort of worked but wasn't too great and the edge it would pick out was 2-3 pixels wide.
I then just tried all of the hydrology tools (ArcMap) and slope as I thought they might help pick out these sharp contacts.  That worked like the above filtering but gave better results when I used the 'Drop' raster from the 'Flow Direction' tool. See below:

This isn't perfect but I thought I could use this, convert to polylines, manually digitise the missed parts, convert to polygons and then clip each sub-set.  I'm not even sure that this would work but converting to polylines gives me the below issue.  Using the simplify lines check box just made the geoprocessing run for almost 24 hours before I cancelled it.

Any ideas where to go from here?


texts: TypeError: list indices must be integers, not str when trying to calculate new records in several fields

texts:
I am trying to calculate new records in several fields by looping through each field in a list and dividing the value of each record by its shape area. My code is as follows:

inFeature = outFeature
my_list = ['a', 'b', 'c']
divfield = "Shape_Area"

with arcpy.da.UpdateCursor(inFeature, my_list) as cursor:
    for row in cursor:
        for n in range(len(my_list)):
            row[n] = row[n] / row[divfield]
        cursor.updateRow(row)

I am getting the error: TypeError: list indices must be integers, not str

All of my fields are stored as double.


texts: QGIS segmentation fault related to numpy on Ubuntu

texts:
I'm running QGIS 2.14 on Ubuntu 14.04 and everything's been working just fine until I upgraded a few Python packages via pip this afternoon. Since doing that, QGIS crashes with a segmentation fault when I try to start it. I don't get any error message, but when I look at the system log it says:

  qgis.bin[9382]: segfault at 7efeea8031a0 ip 00007efeea8031a0 sp 00007ffc54ffdfa8 error 15 in multiarray.so[7efeea7fb000+f000]

Numpy was among the packages I upgraded, and I believe multiarray.so is part of numpy. I've tried uninstalling sudo pip uninstall numpy. That leaves the apt python-numpy package (version 1.8). Then when I try to start QGIS, I get an error saying my version of numpy is too old. I've also tried sudo pip install 'numpy==1.10.4' to downgrade, but I still get that same multiarray.so seg fault.

Please help me if you can. I've got a massive deadline looming, and this is a really bad time to be having these sorts of issues.


texts: Updating joined excel file in ArcGIS Desktop when adding new sheet?

texts:
When adding a new sheet to a joined Excel file in ArcGIS Desktop, the new sheet will not show till I exit both arcgis and excel.

Is there a better way to update the connection?

I have some projects based on an excel workbook joined to a specific feature class (based on an identical field). That workbook has many worksheets, each one joined to a different field in the feature class, and sometimes there is the need to add new data, in a new worksheet (I wouldn't want to create a new workbook for each new table that needs to be joined), and create a join. That's where my problem is: I can't find the way to refresh ArcMap, so it recognizes a newly added worksheet in the excel, without restarting ArcMap.


texts: How do i get x, y coordinate values when click on map?

texts:
I am very preliminary in arc gis and in fact i am a web programmer.

I created a local serviece and use this to create a web applcation in ArcGIS > ArcGIS Server for the Microsoft .NET Framework > ArcGIS Server Manager.exe .

Then i copy website from "c:/inetpup/wwwroot/projectName" and open this to VisualStudio.

Now i want when i click on map, display an pushpin on click location and send x, y to server.


texts: QGIS - Merge vector layers - Calculate minimum of common attribute

texts:
I'm a QGIS beginner. I have twenty layers of drive time polygons. There are multiple overlapping polygons between the layers. All attribute tables have identical field names. I can merge these multiple layers without issue, using any of the available tools in QGIS. I have not found a way to retain the minimum drive time for each merged polygon. All of the tools seem to sum by default. Although I have not yet had success rasterizing these layers, my guess is that I can do these calculations with rasters. However, I was hoping there might be a way to do this while keeping everything in vectors. Any ideas?

Example of overlapping polygons between layers:

Layer 1 - Drive time = 160s
Layer 2 - Drive time = 100s
Layer 3 = Drive time = 135s

I want the new merged polygon to have a value of 100s.

[EDIT]

The answer below worked for me but I needed to use the "contains" option, not "intersect" to get the exact results I wanted. Also, I got odd looking layers from my first experiments, with many polygons receiving no values. I traced this problem to invalid polygons, which I corrected with the v.clean function.


texts: Getting an error when using arcpy.management.XYTableToPoint()

texts:
I get this error when I run my script:

Runtime error
Traceback (most recent call last):
  File "", line 15, in
AttributeError: 'module' object has no attribute 'XYTableToPoint'

Script:

import arcpy

arcpy.env.overwriteOutput = True
from arcpy import env
env.workspace = r'U:\skelley\Python\Outputs'
inTable = r'U:\skelley\Python\FWD_GPR_import\FWD_Halsey_WGS1984.csv'
arcpy.management.XYTableToPoint(inTable, "quasars", "lat", "long",
arcpy.SpatialReference(4326))

Can anybody tell me why?


texts: Extract by mask with grid shapefile in ArcGIS Pro

texts: I have a tiff file. I want to extract it by a grid but every shp cell should clip the tiff file to a new tiff file. Perhaps the attached image could make sense.


texts: from/to ellipsode to/from cartesian conversion

texts:
i am trying to convert from/to ellipsoidal coordinates (lat,lng,alt) to/from cartesian(x,y,z)
i wrote the below code to do the conversion, but given the initial values of the gpsLat and gpsLng
i receive certain value ofr the northing and Easting, but the problem is when I convert from
Northing and Easting to get the Longitude i do not get the same longitude which is 8.427421853866282

the equation i am using i found it in

http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion

please let me know how to get the correct values

code

public static void main(String[] args) {
    double semiMajorAxis = 149598023; //in km
    double eccentricity = 0.0167086;
    double radOfCurvature;
    double Northing;
    double Easting;
    double computedLng;

    double gpsLat = 49.01511835412477;
    double gpsLng = 8.427421853866282;
    double gpsAlt = 137.56704112819037;

    radOfCurvature = semiMajorAxis/Math.sqrt(1 - ( (Math.pow(eccentricity,2) * Math.pow(Math.sin(gpsLat), 2)) ) );
    Northing = (radOfCurvature + gpsAlt) * Math.cos(gpsLat) * Math.sin(gpsLng);
    Easting = (radOfCurvature + gpsAlt) * Math.cos(gpsLat) * Math.cos(gpsLng);
    computedLng = Math.atan(Northing/Easting);

    System.out.println("Northing: " + Northing);
    System.out.println("Easting: " + Easting);
    System.out.println("computedLng: " + computedLng);
}

result of the code

Northing: 3.9587600919687316E7
Easting: -2.5566975079396147E7
computedLng: -0.9973561069030983


texts: How do I get a route between two points using PgRouting with TIGER2010 Data

texts:
A little background:

I've got PostGIS 2.0 and all the TIGER2010 data loaded for a U.S. State California.

I've got a table with points stored as data type geometry.  I've learned enough to get the distance "as the crow flies" between two of the points and my samples are satisfactory when comparing them with Google Earth.

Now, I'd like to generate distance via a street/highway route.  The requirements on the route are that it be the shortest.  I only need shortest travel distance as a numeric value.  No pretty maps.

-I've got every TIGER 2010 shape files loaded for California in PostGIS 2.0.  I'm such a newbie it's not clear to me if this is sufficient data for calculating a route using PgRouting.  Is it?

-If TIGER 2010 data is sufficient, does it need to be prepared somehow for a shortest_distance query?

If there's a useful how-to somewhere, then please send that along as a reply.

As a total newbie, I didn't know that the roads data is already loaded in the 'edges' table. Thanks to the commenter below for pointing it out!

Now, follow the instructions here: http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/

That should satisfy the 'lines layer' requirement Fabien is referring to.


texts: How to limit legend to show only the features currently visible in the view?

texts:
Using ArcGIS Server JavaScript API v4.x and the Legend Widget, is it possible to limit the legend to only those features which are currently visible?
For example, in the Legend Widget Sample, only show the green High School classes when zoomed in to this extent:


texts: PostGIS ST_LineInterpolatePoint inconsistent results

texts:
I have a set of items that are a known length along a LineString, so I'm generating a Point and saving in the database using (all geometries are using 4326):

point = ST_LineInterpolatePoint(line.geometry, item.metres / line.length)

I'm then using Leaflet to display the same LineString and a plugin to display distance markers every kilometre. This works: the displayed distance markers match the known length of the LineString.

The problem occurs when I plot the calculated Points of the items. The results are inconsistent and the position of the items are incorrect, relative to the distance markers (not by a consistent error margin).

Can anyone think what would cause this?


texts: QGIS: How to add picture linked to feature in a report layout?

texts:
I'm trying to make a report in QGIS, in the attribute table of my layer, I have a column named "picture1" filled with the relative path of a picture, that illustrate the feature.

In the report, I did a "field group section" with my layer and a map controlled by report. Now I would like to add the picture related to each feature, in order to have one page by feature with its location and its picture. With the button next to "image source", I tried to choose the "field type" that identify each feature and also to write concat(@project_folder,'/' ,"picture1.jpg") in the expression, but neither of them worked...

Do you have an idea how to resolve this?


texts: Defining custom basemap gallery in ArcGIS API for JavaScript?

texts:
I'm using arcgis js spi for 3.13 edition, when I using basemaptoggle widget from https://developers.arcgis.com/javascript/jssamples/widget_basemap.html is ok, but when defining my own basemap is doesn't work, the codes are following below, hopefully some one who can help to solve it.

var basemaps[];
//one step define BasemapLayer
        var basemaps = [];
        var fx_image =new esri.dijit.BasemapLayer({
            url: "----"
        });
//second step define BaseMap
            var imageBaseMap1 = new esri.dijit.BaseMap({
                layers: [fx_image],
                title: "fx_image",
                thumbnailUrl: "esri/images/basemap/terrain.jpg"
            });
            basemaps.push(imageBaseMap1);
//another basemap
            var fx_sl = new esri.dijit.BasemapLayer({
                url: "----"
            });
            var imageBaseMap2 = new esri.dijit.BaseMap({
                layers: [fx_sl],
                title: "fx_sl",
                thumbnailUrl: "esri/images/basemap/hybrid.jpg"
            });
            basemaps.push(imageBaseMap2);
     //third step define BaseMapGallery
            var basemapGallery = new esri.dijit.BasemapGallery({
                showArcGISBasemaps: false,
                basemaps: basemaps,
                map: map
            }, "basemapGallery");
            basemapGallery.startup();
            basemapGallery.on("error", function (msg) {
                console.log("basemap gallery error:  ", msg);
            });

Sorry for the unreadable code, here is the normal sequence

var basemaps[];
//one step define BasemapLayer

var fx_image =new esri.dijit.BasemapLayer({url:"http://services.arcgisonline.com/ArcGIS/rest/services/Specialty/DeLorme_World_Base_Map/MapServer" });

 //second step define BaseMap

var imageBaseMap1 = new esri.dijit.BaseMap({layers: [fx_image], title:"fx_image", thumbnailUrl: "esri/images/basemap/terrain.jpg"});

basemaps.push(imageBaseMap1);

//another basemap

var fx_sl = new esri.dijit.BasemapLayer({url:"http://services.arcgisonline.com/arcgis/rest/services/Specialty/Soil_Survey_Map/MapServer"});
var imageBaseMap2 = new esri.dijit.BaseMap({layers: [fx_sl],title:"fx_sl",thumbnailUrl: "esri/images/basemap/hybrid.jpg"});

basemaps.push(imageBaseMap2);

//third step define BaseMapGallery

var basemapGallery = new esri.dijit.BasemapGallery({showArcGISBasemaps:false,basemaps: basemaps,map:map }, "basemapGallery"); basemapGallery.startup();

basemapGallery.on("error", function (msg){ console.log("basemap gallery error: ", msg); });


texts: Extract from a linestring layer, groups of linestrings to a multilinestring objects layer

texts:
I've got a linestring layer with thousands of linestrings. I want to extract groups of linestrings to individual multilinestring objects.

In the following picture you can see a part of the layer with a lot of linestrings. Each object inside the red circles are composed of many single linestrings (different records in the layer). I want to obtain a layer with one unique multilinestring record for each object in the red circle.

I assume that i will have to control+select the linestrings that compose one object and do something to convert it to a multinestring record and add it to the new multilinestring layer but I can't see how to do this in QGIS.

I am open to using different software if there is not a solution to this in QGIS. I also use AutoCAD Map.


texts:
I am fairly new to GIS in general.

I am having some issues converting a GeoTIFF into an numpy array. There are multiple questions online about how to do it, I've tried, many, with the same result each time.

I have a GeoTIFF showing the repartition of a parasite in a country, which, when I open it (with rasterio.open) shows:

Basically what I'd like to do, is turning this heatmap into multiple polygons by different thresholds, which would turn into some sort of topological map. I would get polygons where the concentration of the parasite would be >0.8 for example, then where it is 0.8>concentration>0.5, etc. I thought the best way would be to convert the image into an array, and simply using np.where(), keep parts of the images with pixel values in between 2 thresholds.

However, when I try to do so, this is what I get:

Numpy turns everything outside of the boundaries of the GeoTIFF into -9999, which I don't mind. However, I lost all the information inside the boundaries, and I don't really know how to keep it. Here is the code used:

data = gdal.Open(PR_files[10])
img = data.GetRasterBand(1)
raster = data.ReadAsArray()
show(raster)

(Note that I have multiple images over the years, with PR_files being the list where I keep all the paths to these images, so I just fetch one as an example. My images only have 1 band by the way).

What I would ultimately like to do, is turn these shapes into a shapely file, and use these Polygons to do some LULC prediction. I have similar images for each year for the past 17 years, and I'd basically try to see if it would be possible to predict the 18th with this method. I have no idea if this would work, but that's why I want to give it a try. This is just a personal project I'm doing to mostly get more familiar with GIS in general.


texts: MapInfo 10.5.1 Update and MS Access Database

texts:
I am having problems with a database that I have developed. The data uses MapInfo for any of the spatial parts of it.

I installed the update to MapInfo 10.5.1 and I now get an error when I
am trying to access an ADODB.recordset. object which previously
worked. I get the error:

Method 'Connection' of object '_CurrentProject' failed

when the database tries to run the code:

rs5.ActiveConnection = CurrentProject.Connection

(where rs5 is a ADODB.Recordset):

Has anyone else seen anything like this? Is this related to the
updated drivers for Office 2010 that were installed with the update? I
am using Access 2007. The database still works on machines that are
10.5.


texts: GeoServer pregeneralized layers not working as expected

texts:
I've set up GeoServer 2.1.4 and installed the feature-pregeneralized-2.1.4.jar plug in and processed my shape files into each of their separate folders with the generalized shapefiles and I've check these and they are correct.

I've followed the tutorial and set up the shapefiles and configured the XML file as per the article, but I'm not getting the expected log file entries.  I'm getting the fact that the datastores are getting loaded from each of the shapefiles, but when I do a layer preview I don't get the switching between the different shapefiles.

Below is the geninfo_shapefile.xml that I have sitting in the root with the folders:

<?xml version="1.0" encoding="UTF-8"?>
    <GeneralizationInfos version="1.0">
        <GeneralizationInfo dataSourceName="file:data/Local_Authority/0/LocalAuthority_OSBL_May2012.shp"  featureName="LocalAuthority_OSBL_May2012" baseFeatureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom">
            <Generalization dataSourceName="file:data/Local_Authority/1.0/LocalAuthority_OSBL_May2012.shp"  distance="1" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/2.5/LocalAuthority_OSBL_May2012.shp"  distance="2.5" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/5.0/LocalAuthority_OSBL_May2012.shp"  distance="5" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/25.0/LocalAuthority_OSBL_May2012.shp"  distance="25" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/100.0/LocalAuthority_OSBL_May2012.shp"  distance="100" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/500.0/LocalAuthority_OSBL_May2012.shp"  distance="500" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
            <Generalization dataSourceName="file:data/Local_Authority/1500.0/LocalAuthority_OSBL_May2012.shp"  distance="1500" featureName="LocalAuthority_OSBL_May2012" geomPropertyName="the_geom"/>
        </GeneralizationInfo>
    </GeneralizationInfos>

When I go to map preview and then check the log I only have:

2012-10-17 17:50:33,961 INFO [geoserver.wms] -
Request: getServiceInfo
2012-10-17 17:50:33,961 INFO [geoserver.wms] -
Request: getMap
    blah
    blah
    blah

Can anyone assist me with this?


texts: Arcpy standalone script: Tool is not licensed error

texts:
I'm trying to execute a standalone python script using arcpy and spatial analyst extension. I've arcmap installed with a trial license.

The script is the following:

import arcpy,arcinfo
import sys
import os

def listFcsInGDB(dataset):
    ''' List all Feature Classes in a geodatabase, including inside Feature Datasets '''
    fcs = []
    for fds in arcpy.ListDatasets('','feature') + ['']:
        if (fds == dataset):
            for fc in arcpy.ListFeatureClasses('','',fds):
                fcs.append(os.path.join(fds, fc))
    return fcs

if arcpy.CheckProduct("ArcInfo") == "Available":
   print "available"

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")
    print "License is available"

arcpy.env.workspace = "C:\\Users\\user\\\Desktop\\input.gdb"

fcs = listFcsInGDB("SPlit")

for fc in fcs:
    try:
        logFile = open('Log.txt', 'a+')
        logFile.write("Create TIN...\n")
        # Process: Create TIN
        outputFC = fc+"_TIN"
        arcpy.CreateTin_3d(outputFC,
                       "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]",
                       fc,
                       "DELAUNAY")
        logFile.write("TIN successfully created. Output FC: " + outputFC + "\n")
        logFile.close()

...

When I run the script on my machine, i get the error "ERROR 000824: The tool is not licensed", while running the same script on another machine with a regular license does not return any error.

The script prints "License is available" and running the createTin tool from arcmap does not return any error.

I think the issue is due to the fact that I have old expired trial licenses registered in arcGis administrator, and my standalone script is not able to retrieve the correct (unexpired) license. Does this make sense? How can I make it work?


texts: Can adding of attribute table fields cause geoDB replica fail

texts:
We have a arcsde replica between 2 sql server geoDB, it is one way replica and our side is the child. Recently our scheduler to synchronize the replica failed.

The parent side add fields in some layers attribute table, but they insist only the removal of the attribute fields can cause the replica synchronize fail, additional of fields should not affect the replication.

Before I compare the replica schema, I want to confirm whether the adding of attribute fields can cause the replication fail.

Both the parent and child side using arcsde 10.1.


texts: Changing the distance unit from degrees to meters on my semiviriogram in the geostatistical tool in ArcMap?

texts: How can I change my coordinates from decimal degrees to meters?


texts:
I'm wondering if someone could point me in the right direction to force the BBOX of OpenLayers.Control.GetFeature to be always in long/lat coordinates. I'm trying to get this with OpenLayers.Class, so the new control will inherit OpenLayers.Control.GetFeature with just the small modification in BBOX. This is what I've done:

OpenLayers.Control.myGetFeature = OpenLayers.Class(OpenLayers.Control.GetFeature,{
        selectBox: function(position){
            var opts = OpenLayers.Control.GetFeature(position);
            var baseSRSutm = this.map.getprojectionObject();
            var layerSRSdd = new OpenLayers.Projection('EPSG:4326');
            opts.params.bbox = this.map.getExtent().transform(baseSRSutm,layerSRSdd).toBBOX(null,firstLayer.reverseAxisOrder());
        }
        return opts;
    });

but I'm getting that some "}" are missing. Is this the right way to use OpenLayers.Class to do what I want?

Any help is much appreciated.

PS. I'm using OpenLayers 2.11

EDIT:

Moving "return opts;" inside the function solved the above problem but it complained about CLASS_NAME is missing. Adding that, it still complained about the same, but after a couple of clicks or boxes it sent the request, however the sent bbox is still in UTM, where could the problem be? here it is the code:

var mycontrol = function(){
        var protocol = new OpenLayers.Protocol.HTTP({
            url: 'http://www.....',
            format: new OpenLayers.Format.GeoJSON({
                ignoreExtraDims: true,
                'internalProjection': new OpenLayers.Projection('EPSG:900913'),
                'externalProjection': new OpenLayers.Projection('EPSG:4326')
            })
        });

        OpenLayers.Control.myGetFeature = OpenLayers.Class(OpenLayers.Control.GetFeature,{
            selectBox: function(position){
                var opts = OpenLayers.Control.GetFeature(position);
                var baseSRSutm = this.map.getprojectionObject();
    console.log('this is the baseSRSutm' + baseSRSutm); // doesnt appear in firebug
                var layerSRSdd = new OpenLayers.Projection('EPSG:4326');
                opts.params.bbox = this.map.getExtent().transform(baseSRSutm,layerSRSdd).toBBOX(null,firstLayer.reverseAxisOrder());
    console.log('this is the opts.params.bbox' + opts.params.bbox); // doesnt appear in firebug
                return opts;
            },
            CLASS_NAME: "OpenLayers.Control.myGetFeature"
        });

      return new OpenLayers.Control.myGetFeature({
            protocol: protocol,
            box: true,
            click: true,
            single: false,
            clickTolerance: 10,
            eventListeners:{
                // some stuff
            }
        });
    };

because the console.log doesn't appear in firebug, I think the problem is in the OpenLayers.Control.myGetFeature and OpenLayers.Class. It seems to me that the "selectBox" function is the one I need, but the subclass itself is not read by the return.


texts: Could not load library "C:\...\ rtpostgis-2.2.dll"

texts:
I've installed PostgreSQL with Postgis (tried both version 9.4 and 9.5 ) on a remote desktop running Windows 7 64bit.

However, when I executed this SQL command: "CREATE EXTENSION postgis" (on PgAdmin run as administrator), this message came up:

  ERROR:  could not load library
  "C:/PostgreSQL9_5/lib/rtpostgis-2.2.dll": The specified procedure
  could not be found.

  ********** Error **********

  ERROR: could not load library
  "C:/PostgreSQL9_5/lib/rtpostgis-2.2.dll": The specified procedure
  could not be found. SQL state: XX000

I can confirm that "rtpostgis-2.2.dll" was there.

Please tell me how to fix this. It's so frustrating to uninstall and then reinstall every time but the problem still persisted.

P/s: On my Home PC running Windows 10 64bit, PostgreSQL is running like a charm. Postgis and Pgrouting extension can be created easily.


texts:
How would I determine what the red point (phi', theta') is for a camera oriented at a known (phi, theta) is from a given image?

The point of this is to determine what (dphi, dtheta) is so that I can rotate my camera to that location.

I was thinking of a equirectangular or Mercator projection, but couldn't quite figure out how to convert map image coordinates to (phi', theta').


texts: Calculate the latitude,longitude 1000 kilometers east from a given location

texts:
I am trying to calculate what would the Google maps latitude and longitude be, for a given location, 1000 east from it. I am not looking for an application which can do, this, but rather a manual way, by using equations.

I found out that there is a way of doing this, if one would know the initial bearing Angle:

bearingAngle = 90
d = 1000 kilometers / 6371 kilometers
x = asin( sin lat1 * cos d + cos lat1 * sin d * cos bearingAngle)
y = lon1 + atan2( sin bearingAngle * sin d * cos lat1, cos d - sin lat1 * sin x)

(Equations are taken from: movable-type.co.uk/scripts/latlong.html, section "Destination point given distance and bearing from start point")

However, the upper equations are giving the following result:

As you can see the "x,y" location has a bit lower latitude then starting "lan,lon" (Paris).
Basically this is what I am trying to achieve:

The problem is that, to do this, I need to know the initial "bearing angle", which I do not.

How can I solve this problem?