texts: google earth engine: query the image collection for valid 'scenes'?

texts:
I would like to be able to do computations based on 'scenes' using the style of the Landsat Algorithms example in the EE Developer's Guide:

var raw = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_001004_20140524');
Map.centerObject(raw, 10);
Map.addLayer(raw, {bands: ['B4', 'B3', 'B2'], min: 6000, max: 12000},
'raw');

from https://landsat.usgs.gov/wrs-2-pathrow-latitudelongitude-converter , I find that my nearest scene center is: path=43, row=36

I can get date_range from my collection with: ee.DateRange(), but I realize that my scene may not exist for date[0] or date[1].

I believe I correctly build my scene ID as:

'LANDSAT/LC08/C01/T1_SR/LC08_043036_20180406'

but apparently this image is not in the collection.

Image.load: Image asset 'LANDSAT/LC08/C01/T1_SR/LC08_043036_20180406' not found.

(same result for each day decrementing backwards from the collection Last Modified date (2018-04-27) advertised by the collection info-sheet; I also flipped the path and the row, no help)

var raw = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_001004_20140524');
// works

question: is there a way to query the collection to determine the number of date-stamped images for a specific scene ID?

I don't see something like ee.sceneImageGrabber(collection, path, row) returning a list of dates/images....

or ee.polygonImageGrabber(collection, polygonID)
also returning a list of dates/images would be nice....


texts:
I am trying to convert the feature classes in data.gdb to a new PCS per the code on this webpage:

import arcpy
import os
arcpy.env.workspace = "G:/W135thStCopy/data.gdb"
arcpy.env.overwriteOutput = True
outWorkspace = "G:/W135thStCopy/data_WebMerc.gdb"
try:
    for infc in arcpy.ListFeatureClasses():
        dsc = arcpy.Describe(infc)

        if dsc.spatialReference.Name == "Unknown":
            print ('skipped this fc due to undefined coord sys: ')
            print (infc)

        else:
            #Determine the new outpurt feature class
            outfc = os.path.join(outWorkspace, infc)

            #Set output coord sys
            outCS = arcpy.SpatialReference('WGS 1984 Web Mercator Auxiliary Sphere')
            break
            #Run project tool
            arcpy.Project_management(infc, outfc, outCS)

            #Check messages
            print(arcpy.GetMessages())

except arcpy.ExecuteError:
    print(arcpy.GetMessages(2))

except Exception as ex:
    print(ex.args[0])

ERROR 999999: Error executing function.

The error seems to be at the call to SpatialReference() since I don't get the error when I insert the break above that line. I have tried removing the underscores in the coordinate system name and still get the error. How can I fix this?


texts: No longer able to change shapefile field types in QGIS 3.6 or 3.8 using Refactor tool

texts:
I'm looking to merge 3 shapefiles and I am changing field names and types to coordinate before the merge process.

I have been able to successfully use the refactor tool to change field types before (integer64 to integer). I am trying to change a string field to an integer, but each time I try to change it the refactored file still has that field as a string. I have also tried to change an integer64 to an integer in this process, and that also is not happening anymore. Changing field names in the refactor process also doesn't work.

Troubleshooting so far - rebooted, tried in QGIS 3.8 (file originally edited in 3.6), and reloading the shapefile layer. No settings were changed in 3.6 that I know of before this started occurring.

UPDATE: Current work around is to create a new field of the desired type then copy information into it. Very time-consuming and would much rather use the Refactor tool if possible.

Sample shapefile: https://drive.google.com/open?id=1U7uu7LHToiX5OwATipDqcZjk9bmeQcLF

The example offending field is i_ZIPCODE.


texts: Drawing a polygon with a hole in it using Mapstraction

texts:
I am using Mapstraction to make a web-based mapping system. The documentation has some serious holes in it. I would like to use Mapstraction to draw polygons with holes in them. Serious holes.

Can anybody tell me how I can do this?


texts: Geoprocessing Results window does not open, as well as bottom right pop-up window with report link (eg: ANN .html file)

texts:
Suddenly, totally out of the blue, the Geoprocessing Results window doesn't open anymore (ArcMap 10.8. I click on the drop-down menu and nothing happens.
The pop-up mini window which appears at the end of the processing usually on the bottom right screen does the same thing: I click on it and it doesn't open to the report file (eg: ANN).
For ANN I can retrieve the report into the ArcGIS folder, but I do need to see the contents in the Geoprocessing Results window for many other things.
I haven't changed any setting.
I am not using two screens.
I have tried cascade windows on the taskbar to check if the windows were open and hidden somewhere.
I have checked that the pop-up windows in the browser were not blocked (just in case).
I am not running the tools from a script or a model or an external toolbox.
The version I was using was 10.8 with 64bit background (it was working fine on it until 4 weeks ago). I uninstalled it and reinstalled it, but the problem was still there. I upgraded to 10.8.1 because I was advised to do so, but the problem persists.


texts: Getting layer label from ArcMap Table Of Contents?

texts:
How do I get a layer´s label from the TOC?

I can´t find a similar question and ArcGIS Help tells nothing about it.  I have searched even in other knowledge bases but with no success.

I need to get a layer´s label (not its description) to be able to work with it further.

Getting a layer´s name is easy but getting its label that appears in the TOC not.

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, data_frame=df):
    layer_name = lyr.name
    layer_label = ????     # there is nothing like lyr.label


texts:
I am using QGIS 2.18 on my old laptop, but I just got a new laptop and don't have time at the moment to learn QGIS 3.  Googled downloading 2.18 but the responses seemed to require me to download some other software before downloading 2.18, which I am suspicious of.
Eventually I will dive in to 3, but for right now what is a safe and reliable way for me to download 2.18 on my new laptop?


texts: Selecting Polygons that have certain raster data within them

texts:
I am using a raster layer with either 1 or 0 in each cell. On top of this layer, I have a European countries shapefile. The countries that contain a '1 cell' in the raster beneath are the Countries I am interested in studying.

Is there anyway to select only the country polygons that have a '1 Cell' in them from the raster beneath?

I am using ArcGIS 10.3 for Desktop


texts: How to change the layer selection view from Radio button to Image in Leaflet

texts:
I have implemented the below example to show the multiple base layers of map in my application using LeafletJs.

Leaflet JS LayerGroup example

I'm able to successfully implement the same but I'm trying to show same in a single tile image option.

Like for example I have two layers
1. Street Layer
2. Gray scale
With the above implementation in getting the below options on the map

Instead I'm trying to implement something like the the one below

So that when the user clicks on the small map icon, it should change the view to the other Layer. As in my case there will always be only two layer.

I have tried to add the pane option but that doesn't seem to be the right option for this behavior. Also I have tried some of the properties mentioned in the below link.

Layer group documentation

Is this feasible with the existing components of Leaflet?


texts: Convert a *.mif format file to a *.csv point file without loosing attributes

texts:
I tried to convert from .mif to .csv in two steps.

Step1- Save as .shp from Layer
Step2- Convert Polygon/line vertices to Points.

Now, when I tried to do "Point Statistics for Polygons", then an attribute is missing in the Points, which I need.

Is there a way to convert exactly all the attributes of the .mif apple to apple to point so that I can use the function "Point Statistics for Polygons"?


texts: How to Choose Geoserver and GeowebCache Configuration Options?

texts:
I have two Win 2008 R2 Servers currently dedicated to running one instance of Geoserver each. We are looking at setting up geowebcache.
Win Servers have 2Gb Ram and they do nothing but run Geoserver at the moment.
For practical purposes, I have 1Mb, maybe 1.2, to allocate to TomCat.

First, does anyone have any experience on how much Ram can be allocated to Geoserver/geoWebcache withoput impacting Windows Performance?

Second, which way should I set it up?

Option 1: One server running GeowebCahe and one server running Geoserver.
Option 2: Geoserver and geowebcache on each server.

Third, the Geowebcache documentation recommends allocating at least 256 Mb of heap memory. It would be reasonable to assume the more memory the better. If Option 2 is chosen, how should I apportion the available memory between Geoserver and GeoWebCache?

Fourth, am I better off running Geoserver in TomCat or in the environment it is distributed in?

Last, are there any pros and cons on performance for a JDK vs a JRE environment?


texts:
QGIS 2.18.

Problem: I had use group stat plugin in order to sum colums and obtain a value, which I saved in a csv-file. I then joined the csv-file with my original Layer in order to show in my map the summarized values i had obtained. Naturally, my Layer has more Objects than the summarized table, resulting in the repeatedly listing on the 5 unique values within my layer. I would like to be able to show these unique values in the print composer, its not working though. I'm able to show my attribute table, but i can't tell him to only show the unique values. Do you have a solution?


texts: What type of data is available at the "census block" level?

texts:
I'm curious what is available at the 'Census Block' level,  the finest scale.  I found this dataset:  http://www2.census.gov/geo/tiger/TIGER2010/TABBLOCK/2010/

It contains HOUSING10 and POP10,  but I would like to get Vacant vs. Occupied Housing and income if possible.

In spite of the fact that it says this is answered elsewhere, it is not.   Only # of houses and total population are addressed in the other post.   The question here is...  What addtional attributes are available at the 'census block' level?  How do you obtain these additional attributes?


texts: Offset a Polygon defined in Lat Lng points by a given distance of miles or km

texts:
I am trying apply a positive offset/buffer (enlarge a shape equally from the borders) to a polygon that is made up of a series of latitude/longitude coordinates by a distance that is defined in either miles or kilometers, as apposed to a distance specified as a magnitude of the shape.

My research has dug up a number of libraries that will enlarge a polygon by a given magnitude, however this magnitude is relative to the units that the polygon is specified in which isn't helpful in the case of lat/lng coordinates given the shape of the Earth and the changing proportions depending on position on Earth.

I need to accomplish this server side preferably and have dabbled with arcgis's rest service which isn't working for me. Also have looked into ClipperJS and gotten close, however couldn't translate it's delta variable, the scale of increase into miles or km.

Can anyone tell me how this might be achieved? I imagine this is not a straightforward problem but certainly one that has needed to be solved before.


texts: How to display only one representative contour line in legend?

texts:
I have a layer of topographic contour lines, graduated by color (using QGIS). My legend wants to display every level and its color. This is not feasible. I need to simply show one representative line, which I'll edit to be labeled "Topographic contour lines: 5-meter intervals". How can I include one representative line, but not all of the various levels?
Thanks
TVZ


texts:
What steps to I have to take in ArcMap (ArcEditor) to input a lat / long and return nearest address?

I found the Reverse Geocode tool but I can't figure out what to put in the "Input Features" and "Output Feature Class" fields or where to input my lat / long.


texts: Clip Raster with polygon if different SRID (Postgis)

texts:
I have a raster with a SRID of 4326 and a polygon of SRID 2275. The raster was uploaded to postgis with 90x90 tiles using the following raster2pgsql cmd...

raster2pgsql -s 4326 -I -M -b 1 nlcd_2011_clip_wgs84.tif -F -t 90x90 public.nlcd_2011_clip_wgs84 | psql -h localhost -U postgres -d routing

I need to clip raster to extent of polygon (different SRID).

I am having some success using the following postgis query...

DROP TABLE IF EXISTS nlcd_clip_2275;
SELECT
rid,
rast,
filename,
st_intersects(polygon.geom,
    ST_Transform(raster.rast,2275))AS result
INTO nlcd_clip_2275
FROM
public.nlcd_2011_clip_wgs84 as raster,
public.buffer2275 polygon;
DELETE FROM nlcd_clip_2275 WHERE result = false

But the processing is slooooooowww!

Is there a faster postgis logic for both projecting and clipping the raster?


texts:
I have seen many GNSS/INS sensors come with "RTK GNSS receivers" which enable you to achieve RTK. I have been looking at the sensor I have got which claims to be RTK capable, the user manual says it only needs GNSS antenna to be set in the rover mode, InertialSense.

Hence my question is: how is RTK GNSS receiver different from GNSS receiver? How can my sensor achieve RTK without an RTK GNSS receiver?


texts: Can Identify tool to return the raw values of rasters?

texts:
I have a raster in ArcGIS 10 SP4 set to bilinear resampling in its display properties.  When using the identify tool it returns the resampled value instead of the raw value.

Can we get the identify tool to natively return the raw value?


texts: Show sensor node location in QGIS MapCanvas

texts:
I have a Postgres db storaging data from a mobile sensor node, and want to show the location in a QGIS MapCanvas.
Data is in a column in a jsonb format (Location and temperature).
I'm trying to use a virtual layer to extract the data from Postgres db layer (not using PostGIS, I'm keeping the original db with no geometry data).
There is no way to extract the data from json column manually, and I can't find any solution in PyQGIS. When I create a virtual layer using a query, the layer shows an error.
vlayer = QgsVectorLayer( "?query=SELECT (object ->> 'longitud')::real AS longitud FROM device_up", "vlayer", "virtual" )

object is the jsonb column. Has anybody worked with something similar?


texts: Technology moves on but this is still a relevant Question How to embed a map into Access Database?

texts: I'm currently trying to build a front end interface for a land-parcel management system in Microsoft Access and in displaying individual parcels with their attribute data - I wanted to display on the form a relative map of the parcel on OpenStreetMap - I've come across the MapWinGIS https://www.mapwindow.org/ which I can use to add an active x controller to the form in Access - I'm missing the bit were I can link a geom point, polygon or line from a table (the table is an external table linked through the PostgreSQL ODBC driver from a PostgreSQL PostGIS enabled database...


texts: Remove dark background around clipped image in QGIS 3?

texts: I used QGIS 3 for clipping a satellite image using a shapefile.  The output raster gave a clipped image with dark background. How to fix this issue?


texts: How can I refer to a model variable in a Field Calculator or Pyculator formula?

texts:
I am new to QGIS 2.0.1 Modeler (Sextante) and not very able in Python programming.

My question regard the use of field calculator (or advanced python field calculator) inside the modeler.
Is it possible to refer in the field calculator formula to an input variable of the sextante model?

Example for Buffer:
In the graphical modeler I define 2 model inputs:

layer = the point layer to be buffered (type: vector layer)
buff_dist = the buffer radius (type: number)

I use the "Fixed distance buffer" algorithm, setting "Imput Layer" = *layer* and "Distance" = *buff_dist*. It creates a temporary output shapefile.
I want to add a new numeric field in the output of the "buffer" algorithm, with value equal to the buffer distance previously defined (buff_dist).

How can I refer to the variable buff_dist in the Formula (or Global Expression) of the "Field Calculator" (or "Advanced Pyculator") algorithm??


texts: In QGIS is there a way to make one polygon out of several

texts:
I'm trying to use QGIS to get a more simplified vector country data. For example Canada has several bits to the north all as separate polygons where as I'd be happy to just take a rough outline of the whole of it as one polygon.

I looked at using a convex hull, but that's a little too rough


texts: Delination of half basins or hillslopes from a flow direction map or/and a stream network

texts:
I am using the GrassGIS r.stream.extract extension to delineate a stream network with the options 'mexp' (Montgomery exponent) and 'd8cut' (Use SFD above this threshold).
It works perfect so far and is a great improvement over r.watershed. My problem is that if I use r.stream.basins afterwards I cannot extract half basins (left and right site of a river) as it was possible in r.watershed. I wonder if there is a tool with whom I can do that? r.stream.extract delivers a flow direction map and the stream network.

Best Regards


texts: Exporting NetCDF layer with defined colors using QGIS?

texts:
I have average wind speed NetCDF file (.nc) that I imported to QGIS. Then I defined colors for wind speeds in equal interval (3, 3.5, 4, 4.5 etc). Now I would like to save this colored map for sharing purpose. When exporting and saving as GeoTiff, the colors change to black and white.

What should I do?


texts: Network Analyst does not take oneway street from attribute table

texts:
I want to build a network dataset to calculate travelling times for buses (service area). I have on shapefile with the streets from my home town, which has a column intitled 'ONEWAY', with values of -1 (to-from), 0, 1 (from-to), but the network is created without the oneway restrictions on those streets. Even if i try a route, for example, it is generated without taking in consideration the oneway street.
Can anyone help me understand why or what i could do?


texts: Failing to Map Reducer using Google Earth Engine

texts:
I'm trying to normalize a collection of images, each to it's own mean value (over a certain region), but I keep getting an error inside the map() function that I can't find anywhere else in documentation or online.

Internal error: function argument not initialized.

Here is my code:

var col = ee.ImageCollection('COPERNICUS/S1_GRD')
  .select(['VV'])
  .filterBounds(PolyOfInterest)
  .filterDate(StartDate,EndDate);

var normalize = function(image) {
  var tmp = image.clip(PolyOfInterest);
  var stats = tmp.reduceRegion({
    geometry: roi,
    reducer: ee.Reducer.mean(),
    bestEffort: true,
  });
  print(stats); // I don't know how stats will be structured, so I haven't yet applied it to the image in the next line.

  return image; // Scale the image so that the mean value is mapped to 1 (values above the mean are clipped)
};

var scaledCollection = col.map(normalize);

Has anyone else encountered/solved this problem before, or could help me?

The problem persists even if the entirety of the normalize function is replaced with simply print(image);.


texts: Clicking through one layer onto another in Leaflet?

texts:
I have a map with a polygon layer and a marker layer.  They are both derived from GeoJSONs, which are identical in every way except for the geometry (the marker layer is simply the centroids of the polygon layer).  I currently have some functionality for the user to click on the polygon layer, and the selected polygon is coloured, added to some tracking arrays, and some summary information derived from the GeoJSON layer is displayed elsewhere.

The problem is that the associated marker layer is quite large, and makes it annoying to click on the polygons.  (See image below)

Clicking on the markers does not count as a click on the polygon layer.  Rather than having to write some new functions to deal with grabbing values from the marker layer and looking them up in the polygon layer to do the same thing, I would rather just somehow "click through" the marker layer directly onto the polygon layer underneath, or somehow forward the click coordinates to it.

I am aware of the fireEvent method in Leaflet, but I cannot seem to figure out a way to use it to send the click coordinates to the polygon layer. Has anybody done something like this before, or perhaps have a better approach?

By the way, I have already seen Trigger a click event on a leaflet map? which appears to be the same, but if you read the accepted answer, the author was simply overcomplicating a routine task, and so the question of "triggering a click event with coordinates" was never actually answered

This is the code I am working with now, it doesn't work, but I can't find any working examples of passing a click to a layer:

onEachFeature: function(feature, layer) {

            layer.bindTooltip(feature.properties.name, {
                sticky: true
            });

            layer.on('click', function (e) {

                console.log(e);
                var lat = e.latlng.lat;
                var lng = e.latlng.lng;

                var latlngPoint = new L.LatLng(lat, lng);
                console.log(latlngPoint);

                markerPolygonLayer.fireEvent('click', {
                    latlng: latlngPoint,
                    layerPoint: map.latLngToLayerPoint(latlngPoint),
                    containerPoint: map.latLngToContainerPoint(latlngPoint)
                });

            });

        },


texts: Removing thin rectangles from a shapely polygon

texts:
I am detecting rooftops and output them as shapley polygons in square, U and L shape. Sometimes the output is in U shape with a very thin rectangle on one side ( true shape is L shape) which needs to be cleaned. I have tried using perimeter/area ratio to detect such parts of a polygon and remove them but I am not able to implement it.
An example is shown in figure below.

I am looking for a solution preferably with shapely.

The polygon in the example in GEOJSON can be found here.


texts: Download Black Sea contours (OpenStreetMap)

texts:
I need to download the contours of the Black Sea. Using the program QGIS (plugin QuickOSM), I was able to download the contours of the Sea of Azov. But I could not do this for the Black Sea. Need a script that will give me the Black Sea multipolygon. I get the Sea of Azov using a script (below) no problems.
https://overpass-turbo.eu/

<osm-script output="xml" timeout="75">
    <id-query {{geocodeArea:Sea Of Azov}} into="area_0"/>
    <union>
        <query type="relation">
            <has-kv k="place" v="sea"/>
            <area-query from="area_0"/>
        </query>
    </union>
    <union>
        <item/>
        <recurse type="down"/>
    </union>
    <print mode="body"/>
</osm-script>


texts: modal resampling in GDAL: what happens when two classes are equally represented

texts:
gdal_translate has a parameter to resample the values based on the most frequent class in the output pixel (-r mode). However, this most frequent class does not necessary exist (when two or more of the most frequent classes have exactly the same number of pixels.

Does anyone know what happens to the output pixel value when two classes with the most pixels are exactly the same ? Does it take the first discovered value, the last one, the smallest one, the larger one, a random one or any other rule?  At least it seems to select a value.


texts:
I use QField for tree surveys. I use Geometry generator to create the crown (canopy) spread based on crown radius taken at N, S E and W. This works, but the crown is 'virtual' and cannot be edited. Rather than use 'Geometry Generator' or a 'Virtual Layer', I would like to automatically create the geometry in an editable layer. I'm aware QField does not yet support Python. But is there some other way of creating editable geometries based on the attribute values of a layer?
O


texts: Is it possible to show 3d models using QGIS?

texts:
I can show vector layers from PostGIS using QGsMapCanvas and QGsVectorLayer in my standalone application. But now I need to show 3d models, for example 3dmax or sketchup, of some buildings over vector layers.

Is possible at all?

If it is possible, how can I do this?

Maybe there are any plugins for this?


texts: Select polygons that contains points in a separate shp file

texts:
I am a very new QGIS user and have years of experience using the basic functions of ARCGIS. I am having trouble using the select by location feature in QGIS. I have a point layer with locations of a species. I also have a polygon layer that is a grid over a particular area. I need to select the grid squares that contain one or more species locations.

I have made sure both the point and polygon layers are in the same CRS. I have tried it with all of the points selected and none of the point selected. The result I keep getting is after a lot of thinking none of the grid squares are selected.


texts: dataset = gdal.open(file) how can I get the compression method?

texts:
I am working with various file sets and one of the things I have not been able to figure out is how to get the compression used in the file that I opened.

The challenge is that I have some different data sources and I want to be able to address the way that I handle the file based upon if the source has compressed the data (lossy) or not.

Does anyone know how to get the compression type from the gdal dataset in python?


texts: Creating path based on waypoints using QGIS?

texts:
I am trying to map flight paths and compare historic paths with current ones (so I later need to perform some calculations as distances between them). For this task I am using QGIS on OSX 10.9.1

As an input, I have the waypoints of every path that I want to connect linear to a line (see picture).

When using the “add delimited text layer” function, I can map the waypoint. However here I am stuck. I already tied the points to path plugin, but I can´t see any line in my map.

I read that it might be possible by transforming the points into a multilinestring in WKT format. Unfortunately, I never worked with WKT files and do not even know how to create them. Also, there is a really long list of waypoints in the actual sample.

All the work should not be done in a “dirty” way (just draw the line), since there are some calculations about distances (e.g. shortest distance between paths, etc) to be performed at the end.


texts:
I am trying to make a custom color ramp by using preset color ramp in ArcGIS 10.5. I need to show 19 colors. The default number of colors that can be input in preset color ramp is 13.

Is there any way to increase the number of this default input value?


texts:
I am using QGIS 3.10 and the integrated GRASS plugin.
I am trying to run r.watershed but I do not understand where I can pass in more parameters & flags to the r.watershed task.
There are tons of them as suggested by the docs (https://grass.osgeo.org/grass78/manuals/r.watershed.html)
See the red box - is this all I can provide to r.watershed via QGIS?


texts:
I'm relatively new to python scripting and I'm trying to batch process a number of files, with the first command being to convert formats.

gdal_translate -of KEA NDVI_20160101.tif NDVI_20160101.KEA

works perfectly when inputting the command directly into the terminal but when I copy and paste the command into a text editor and try and run the script from the destination with:

python ./SCRIPT_NAME.py

I get SyntaxError: invalid syntax


texts: How to generate a parameter defined rotated grid in QGIS, GDAL/OGR and/or Python?

texts:
Following this question which answers partially to my question, I'd like to go further and find out how I could generate a parameter-directed grid in QGIS, with origin X,Y coordinates, Direction/Angle, number of cells, size of cells.

My problem initially came from a need to add grid squares easily to an existing rotated grid. I could do it by hand, but using a script, I thought might later be able to pass on some attributes while creating the additional grid.

All resources I find explain how to create an horizontal grid with square or rectangle like cells, but I didn't find yet how to generate an inclined grid (with an angle and a Origin X/Y), ideally in the form of a QGIS python script with parameters in the processing toolbox for easy sharing.

I think it's just a matter of adaptation of the following script found out by @Gene to include some angle consideration. I'll try to answer my own question when I find the right formulas but if any of you have already some resource I could use or if you would like to take up the challenge... ;)

#Create a memory layer limited here to the geometry of Polygons
class Crea_layer(object):
    def __init__(self,name,type):
        self.type=type
        self.name = name
        self.layer =  QgsVectorLayer(self.type, self.name , "memory")
        self.pr =self.layer.dataProvider()
    def create_poly(self,points):
        self.seg = QgsFeature()
        self.seg.setGeometry(QgsGeometry.fromPolygon([points]))
        self.pr.addFeatures( [self.seg] )
        self.layer.updateExtents()
    @property
    def disp_layer(self):
        QgsMapLayerRegistry.instance().addMapLayers([self.layer])

#Script to create a squared grid layer from an existing layer
from math import ceil
canvas= qgis.utils.iface.mapCanvas()
# first layer
layer = canvas.layer(0)
xmin,ymin,xmax,ymax = layer.extent().toRectF().getCoords()
gridWidth = 1000
gridHeight = 1000
rows = ceil((ymax-ymin)/gridHeight)
cols = ceil((xmax-xmin)/gridWidth)
ringXleftOrigin = xmin
ringXrightOrigin = xmin + gridWidth
ringYtopOrigin = ymax
ringYbottomOrigin = ymax-gridHeight
pol = Crea_layer("grid", "Polygon")
for i in range(int(cols)):
    # reset envelope for rows
    ringYtop = ringYtopOrigin
    ringYbottom =ringYbottomOrigin
    for j in range(int(rows)):
        poly = [QgsPoint(ringXleftOrigin, ringYtop),QgsPoint(ringXrightOrigin, ringYtop),QgsPoint(ringXrightOrigin, ringYbottom),QgsPoint(ringXleftOrigin, ringYbottom),QgsPoint(ringXleftOrigin, ringYtop)]
        pol.create_poly(poly)
        ringYtop = ringYtop - gridHeight
        ringYbottom = ringYbottom - gridHeight
    ringXleftOrigin = ringXleftOrigin + gridWidth
    ringXrightOrigin = ringXrightOrigin + gridWidth

pol.disp_layer


texts: Iterating feature class output file names using ModelBuilder?

texts:
I am trying to iterate a series of feature classes (polygons) in model builder to convert a series of polygons to rasters, however, I am having an issue with the output file names.

After inserting the "iterate feature class" tool into the model builder window, and linking my gdb storing all the polygons as the input, the output feature (green bubble) automatically becomes the name of my first polygon. As a result the model that I attempt to run seems to recur to that named polygon rather than then going to the next polygon in the list. Output rasters are created, but overwrite that file name rather than generating a new raster with a new name that matches the subsequent polygons.

What am I doing wrong?


texts: Why are some tools on Tools toolbar not working in ArcMap?

texts:
I am using ArcGIS 10.2.2 for Desktop, i've managed my map in the layout view, everything done, i was fine with the outcome, exported into a pdf file, then switched back to data view.

Now i am completely unable to use some of my tools from the "Tools" toolbar on this mxd.

If i open another mxd or add a new dataframe on those, everything works just fine, but on this one i am unable to zoom in/out or move it, however i am able to select features.

Also reference scale changed for (seemengly) no reason. I guess somehow i fixed even the dataframe view, not just the layout.

What happened? How can i move my map in data view again? Will the solution have effects in the layout?


texts: Trying to create different sets of random points through ArcPy yet several of the files are duplicates?

texts:
I'm trying to create 50 sets of random point shapefiles in 37 different folders for a total of 1,850 random point shapefiles.  However, upon running the script only a handful of sets are actually random.  For example, out of 50 sets the first 9 would all be identical, then the next 8 would all be identical, and the next 9 would be identical etc.  The function successfully creates a random point file but for some reason makes the exact same file for the next several iterations.  Below is my Arcpy script.

for x in range(4,41):

...     for y in range(1,51):
...   arcpy.CreateRandomPoints_management(filename,str(x)+"kmran"+str(y),"polygonboundary","",157,"","POINT","")

I'm unaware if this is a known bug in either the random point generator tool or arcpy.  I'm using ArcGIS 10.3.1.

EDIT:
Just to be clear this happens regardless of how many runs or points I generate.

EDIT 2:

SOLVED THE PROBLEM!  Turns out that I wasn't giving ArcGIS's random point generator algorithm enough time to create a new set of shapefiles.  I threw in a time.sleep(2) at the end of the 2nd loop to add a 2 second delay beteween each file generation and it worked and generated all unique files!


texts: how to run multiple instances of tile server on a single machine?

texts:
I installed OSM tile server on Ubuntu 14.04.
I followed instruction given here https://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/ .
I want to perform "vertical load balancing" for case of tile server.

In case of Geoserver , it is done like this http://boundlessgeo.com/2013/04/geoserver-in-a-clustered-configuration-part-2/
These steps can't be followed in case of tilserver because tile server consists of many individual components like Mod_tile, renderd, mapnik, osm2pgsql and a postgresql. Mod_tile is an apache module

For running many instances of tile server,I am thinking of running multiple instances of Apache web-server(in which renderd and mod_tile are plugged on which tile request will be received)with HAProxy.
But I can't find exact solution.

Is there any other solution possible?

If running two instances Apache2 is only solution then how to do it?

(I have gone through some sites for many instances of apache2 but necessary httpd.conf of Apache2 is not available in Ubuntu ,so I cannot configured necessary module )


texts: Convert this CSV to a qgis usable format

texts:
I have been given a CSV (which has been converted from a RDF). I want to convert the CSV file to a usable format so that I can display these polygons in QGIS.

The original file and some other formats are available here, but no format seems to work with QGIS. http://datos.santander.es/dataset/?id=distritos

How can I convert it to a usable format, such as GeoJson?


texts:
can somebody help me understand (probably giving visual examples) what does this sentence mean which is quoted from the wiki page of nine intersection model:

  Those points that are in the boundaries of an odd number of its elements (curves).

https://en.wikipedia.org/wiki/DE-9IM

I run the following example in R:

elem1 <- matrix(c(0,0,0,1),2)
elem2 <- matrix(c(0,0,0,-1),2)
elem3 <- matrix(c(0,-1,0,0),2)
elem4 <- matrix(c(0,0,-1,-2),2)

mls <- sf::st_multilinestring(list(elem1,elem2,elem3,elem4))
plot(mls,type='o')

sf::st_boundary(mls)

The output is all 4 points on the boundaries of each individual element.

MULTIPOINT(-1 0, 0 -2, 0 0, 0 1)


texts: Counting Points in Polygon to analyze data using QGIS

texts:
I am using QGIS, and I have a shapefile that has 23 different land classifications. I have another file that lists coyote observations as points.  I would like to find out how many coyote observations are found within each of the 23 categories.
I can create a separate shapefile for each land type and then run "count points in polygon" for each one but does QGIS provide an alternative approach?


texts:
I am using Esri ArcMap 10.3:

I have a very large table so I do no want to do this manually it will take me months (50k+ rows of data)

My table looks something like below example.

Purpose           new field
multi family
apple
jack
random
random

Of course my real table is way bigger and the 'purpose' field has many more unqiue values.

So what I would like to do is use the field calculator on my "new field"
and have it add a value of 'xx' for each row that had a purpose value of 'apple' and jack'.

So basically if purpose was 'jack' or 'random' then new field will be 'xx'

This could be done with python im sure but I was just going to use the calculate field tool and write in an expression.


texts: US Tiger or Open Street Map Route Data: One Way Street & Intersection Data

texts:
We're building a U.S. street routing application (finding paths for cars between two or more locations using public roads) and would like to use open source data US Tiger Data or Open Street Map or other if available nationally.

We have worked extensively with the US Census Tiger Street Segment data before (we built a geocoder http://maplarge.com/geocoder using Tiger), and using the Tiger data is our preference. However, I'm missing two pieces of important data that I think should be there but I can't find them in the documentation:

Question 1: Direction: I know some streets segments are One-way, meaning that traffic is only allowed to flow in a particular direction, and I need to know if the Tiger or OSM have direction attribute data that can be linked to street segments. Specific links to documentation / data / examples would be awesome.

Question 2: Connectivity. When two roads cross it implies an intersection.  However, you can't always go from one road to another.. for example, when a local road crosses a limited access highway, connectivity is often limited to certain "on ramps". I would like to know if the Census/OSM data can reliably provide information about connectivity at intersections. Specific links to documentation / data / examples would be awesome.

Thanks!

(ps this is a high volume application that requires our specific architecture, please don't dodge the question and suggest ArcGIS, MapPoint or a web service.)


texts: How can i use Multiple filter in wfst leaflet?

texts:
I'm using wfst plugin in leaflet. I can send request just by one filter and get response correctly like this:

var layer = new L.WFS({
            url: 'http://localhost:8080/geoserver/ows',
            typeNS: sf,
            typeName: road,
            crs: L.CRS.EPSG4326,
            geometryField: 'geom',
            filter: new L.Filter.Within(values[0], values[1], values[2])},
            style: {
                color: 'blue',
                weight: 2
            }
        }).addTo(map)

How can I send two or more filter Simultaneously like this:
filter: new L.Filter.Within(values[0], values[1], values[2])}
filter: new L.Filter.EQ(values[0], values[1], values[2])}


texts: Given coordinates, unclear how to convert to GPS

texts:
I believe this is "Township and Range" coordinates.
They were handed to me on paper as: 19T 0345856 over 5000288.

I want to plug into Google Earth or my GPS and find.
Apparently they are the location in the middle of Maine, USA to a old cavern. The source is an old man, in his 90's and I could not ask technical questions.


texts: How to merge shapefiles following a particular order?

texts:
I have a dataset containing 171 shapefiles, their names as "shp0_0" to "shp0_171". I need to merge them into a single shapefile. When I use the merge tool and load all the .shp from the folder, it gives me random order, they are not listed by the filename order.

I need the output attribute table's FIDs to correspond with the input shapefiles: shp1 to be in the 1st row and have FID 0, shp2 to be in the 2nd row and have FID 1 and so on.

Is there a way to list all shapefiles to be merged numerically sorted by filename?

I work in arcmap 10.3 (using arcpy).


texts: More creative or efficient solution for a proximity table (1:m) to all points within X distance, keeping obs. where no points fall within X distance?

texts:
I have a proximity question. I am open to Arc or Python solutions. I am mainly asking because I think I am lacking a little creativity in approaching the situation.

The Situation

I have two point files (FileA, n = 100, and FileB, n = 10). I would like to know the nearest distance from each n of FileA to the n of FileB so long as the n of FileA is within X distance of FileB.

Current Approach

I am currently attempting an Arc solution via a 1:m spatial join, and then generating near distances based on the joined object IDs.

The Question/Problem

I'm trying to push my GIS skills and am at a creative roadblock as to how I could improve my solution. I see potential as I may be asking a roundabout buffer question, or perhaps missing some kind of crucial tool/intuition about the problem.

Apologies if this is not really a Stack appropriate question/format - will remove ASAP if so!

Solution via PostGIS (not quite comfortable with PostGIS yet): How to get distance between one point and all other points in a table?


texts: Adjust PictureMarkerSymbol size to distance (ArcGIS 10.2.5 SDK Runtime for .Net)

texts:
I'm working with ArcGis 10.2.5 Runtime SDK for .Net.

I have created a PictureMarkerSymbol that shows data of an area of interest in miles, and the question is, how can i fit the symbol size to a given distance?

I tried with GeodesicBuffer with point geometry, but instead of resize the symbol, change is applied to each pixel of the image, keeping same size and blurring the source image.

I was thinking on an alternative, maybe use polar calcs to get coordinates, measure distance in pixels between two points and resize image with that information, but this is very inefficient and i can't believe there is any better way to do this.


texts:
Is it possible to turn off the 'hover' feature in CartoDB?

I have toggled off all the fields from the 'hover' menu, which in the visualization wizard shows as being 'off', but when I visit the embedded map on the webpage consuming it, the hover is still on, and 'no data available' shows in the hover window.


texts: I want to convert my text file to grib file. But it seems not possible. How can I code grib files in Fortran or C#.


texts: How to create a buffer around a point that takes terrain elevation in account?

texts:
In QGIS, I am creating a few different buffer around points to simulate walking distances. For example, 10 minutes walk about 800m at 5km/h. Using population density, I then clip a population layer which gives me the population that is within a 10 min walk from those points.
I am working on a hilly city where elevation does matter.

How can I adjust the buffer with elevation?


texts: Are there existing maps of U.S. county boundaries?

texts:
I'm relatively new to Carto and very new on the latest interface.

Is there a way to find and join an existing map of Arizona county boundaries to my own data file for the same counties?


texts:
Possible Duplicate:
  Raster incorrectly reprojected to OSGB(27700)

The crs is clearly off when the gridlines run straight and the map lines run at an angle.

I have set the crs to match the project crs, enabled on the fly projections and nothing.

Only thing that happens is it automatically produces a custom crs to produce this off centre projection.


texts: ArcGIS: Compare 'similar' values between multiple tables

texts:
I'm looking to see if there's some code out there I can manipulate, or if I can utilize a model to essentially compare field values between 3 or more tables.  The two issues here are

I have to relate through an intermediate table to go from the Main table to the various tables I want to compare (more below), and
the string values in the main table aren't exactly the same as the values in the compare tables.

Example:
Main_Table has a field Restricted with possible values:

Unrestricted
Restrict 1
Restrict 2

Primary Key field Main ID relates to a Link_Table with fields: Main_ID,  Table_A_ID and Table_B_ID

Table_A_ID and Table_B_ID may or may not be null.  If they have values, then they relate to their corresponding compare tables (Table A and Table B).

Table A and B both have Restricted fields.  But the values here (strings) are either "X" (meaning restricted) or NULL (meaning unrestricted).

Essentially, I want to compare the Restricted value in Main_Table to the Restricted values of Table A and B.  If there are any differences, I want to either select them in Main_Table or output the results in a new table (whatever is easier.)

So, if Main_Table is 'Restrict 1' (or 2), and either Table A or B is NULL, then they fail and that record is selected in Main_Table (or outputted to a new table).

There's quite a bit involved here and I'm not sure if there's already some code out there I can start with.

Update 1 I think screenshots would be confusing, as I'm using generic field and table names here in place of my real data.  But I can give you a snapshot of what the tables might look like.  I should note...record 101 has a related value in Table_B but not Table_A.  I still want to compare Main_Table to Table_B in that case, even though there is no related record in Table_A.


texts: Best way to credit MapTiler and OpenStreetMap in layout manager in QGIS

texts:
I want to use the MapTiler vector tiles plugin and their awesome maps as basemap for some maps. I'm no super experienced doing this and that's why I wondered what the best way is to credit them when moving to the layout view. In the normal map-view in the lover left corner appears this automatically:

However, when I switch to the layout view it isn't there anymore. I could create a textfield and add them manually, but maybe there is an even better way.


texts: ArcObjects: How to create a composite Geotransformation and use with the GeoProcessor

texts:
Using arcobjects, I am trying to programmatically create a composite geotransformation and then use to project features with the  geoprocessor.  Found old VB6 code out there on the web, but I can't get it to work in C# .NET. When I try to use the projection tool in the geoprocessor, it fails and sends back a message saying "Invalid geographic transformation."  The code I am using is below. Either I have built the composite transformation incorrectly, or, the geoprocessor cannot find it.  Any ideas?

EDIT UPDATE:

private void createAndUseCompositeGeoTrans()
    {
        Geoprocessor gp = new Geoprocessor();
        gp.OverwriteOutput = true;
        object sev = 2;
        object dt = "";

        try
        {
            ISpatialReferenceFactory2 spatialReferenceFactory2 = new SpatialReferenceEnvironmentClass();
            ISpatialReference in_spatialReference1 =
                spatialReferenceFactory2.CreateSpatialReference(esriSRProjCSType.esriSRProjCS_NAD1983SPCS_WASouthFT.GetHashCode());

            ISpatialReference out_spatialReference2
               = spatialReferenceFactory2.CreateSpatialReference(esriSRProjCS2Type.esriSRProjCS_WGS1984SphereWebMercator.GetHashCode());
            ICompositeGeoTransformation compGeoTrs = (ICompositeGeoTransformation)new CompositeGeoTransformation();

            //NAD1983_To_WGS1984_1
            IGeoTransformation geoTransformation_A = spatialReferenceFactory2.CreateGeoTransformation((int)esriSRGeoTransformationType.esriSRGeoTransformation_NAD1983_To_WGS1984_1) as IGeoTransformation;
            compGeoTrs.Add(esriTransformDirection.esriTransformForward, geoTransformation_A);

            // WGS_1984_Major_Auxiliary_Sphere_To_WGS_1984
            IGeoTransformation geoTransformation_B = spatialReferenceFactory2.CreateGeoTransformation(108100) as IGeoTransformation;
            compGeoTrs.Add(esriTransformDirection.esriTransformReverse, geoTransformation_B);

            compGeoTrs.Name = "NAD1983_To_WGS_1984_Major_Auxiliary_Sphere";

            IMapGeographicTransformations mapGeotrans = (IMapGeographicTransformations)map;
            IGeoTransformationOperationSet geoTransOperationSet = mapGeotrans.GeographicTransformations;

            geoTransOperationSet.Set(esriTransformDirection.esriTransformForward, compGeoTrs);
            geoTransOperationSet.Set(esriTransformDirection.esriTransformReverse, compGeoTrs);
            ESRI.ArcGIS.Carto.IActiveView activeView = (ESRI.ArcGIS.Carto.IActiveView)map;
            activeView.Refresh();

            ESRI.ArcGIS.DataManagementTools.Project projectTool = new ESRI.ArcGIS.DataManagementTools.Project();
            projectTool.in_dataset = "C:\\projects\\clients\\Freeport\\CAD_Tool\\sampledata\\Point.shp";
            projectTool.out_dataset = "C:\\projects\\clients\\Freeport\\CAD_Tool\\sampledata\\Test_002.shp"; ;
            projectTool.in_coor_system = spatialReferenceFactory2.CreateProjectedCoordinateSystem(in_spatialReference1.FactoryCode);
            projectTool.out_coor_system = spatialReferenceFactory2.CreateProjectedCoordinateSystem(out_spatialReference2.FactoryCode);
            projectTool.transform_method = "NAD1983_To_WGS_1984_Major_Auxiliary_Sphere";
            gp.Execute(projectTool, null);

        }
        catch (System.Runtime.InteropServices.COMException CE)
        {
            //throw CE;
            printMsg(rTxtBx_console, CE.Message + "\n\nLast GP Message:\n" + gp.GetMessages(ref sev), MsgTypes.err);
        }
        catch (Exception E)
        {
            printMsg(rTxtBx_console, E.Message + "\n\nLast GP Message:\n" + gp.GetMessages(ref sev), MsgTypes.err);
        }

    }


texts: Issue with recently cached ArcGIS services

texts:
ArcGIS for Server 10.0 SP5

Virtual Windows Server 2008 R2 Enterprise

16.0 GB RAM, 4 Processors (2.00 GHz each)

I am having issues with two services that I have cached as basemaps for my web applications.

Here's a little back story - These services have been in use for a couple years, and every so often (about 6 months) these caches are rebuilt to address changes in basemap items (i.e. street centerline changes). We use a variety of different basemaps in our web map applications, and this weekend I recently updated the cache on two of them since it takes so long to do them in our current setup (one machine environment).

These services worked just fine, then after recaching them they appear to be fine, until a user spends a little more than 2 minutes on an application using the service. Then the basemap tiles begin to lag and not load. Network debugging shows the pending processes, and SOMETIMES the tiles do load after about 1-2 minutes. MOST times they do not load for at least 4 minutes. Remaining times they just stop loading, and process stay pending until they timeout with a 500 server error.

I have an open ticket with Esri, but they have yet to pinpoint the issue. Some changes we have tried include:

Editing the impersonation for the rest and services config files from true to false
Running post installation
Changing services to low isolation and increasing instances to 8

Results for each change are as follows respectively:

Caused Bad 404 errors for REST page.
This was done in to try and fix above 404 error, it did not. Impersonate values changed back to true to fix 404 error.
No changes seen with services. Same lag issues as before.

In addition to the above, the ArcGISServiceAppPool process periodically will cause the CPU on the server to peg out at 100%. I have seen the KB articles (w3wp.exe, although it is for 9.x versions, and lsass.exe) related to this issue. Implementing the suggested fixes are what caused the 404 issue.

I read through this post and it relates to some of the issues I am also experiencing.  I have checked the log and do not see anything that helps identify the issue at hand.

This question is also posted to ArcGIS forum.  User asked if service caches were being created on demand.  Answer: They are not, however both with and without this option have been tested and the results are similar.

I believe I may have found why the issue is occurring only at some levels.  It appears that some tiles are missing at these levels.  Since some tiles are missing, users still pan/zoom which then become network processes that never complete.  If enough of these add up, it bogs the server down (since most other requests are fulfilled, it's hard to monitor performance).

I think from here I am going to recache the map at the last three levels, then copy the updated levels to the original cache and see if this fixes the issue.


texts: Querying County Open data endpoint and filtering by radius

texts:
I'm using python's requests and arcgis libraries for to query an open data endpoint for my county GIS data.   I'm trying to query the recent sales data within a 2 mile radius of an address.  I've worked out the sql "where" logic to get the "recent" results but can't seem to find a way to do the radius filter.  Here is what I've tried:

from arcgis.gis import GIS
from arcgis.geocoding import geocode, reverse_geocode
from arcgis.geometry import Point
from arcgis.features import *

gis=GIS()
address = "2526 westmoreland dr raleigh nc"
wake_url = 'https://maps.wakegov.com/arcgis/rest/services/Property/Parcels/MapServer/0/'

address_geo = geocode(address = address, out_sr=102100)[0]
address_geo

  {'address': 'Westmoreland Dr, Raleigh, North Carolina, 27612',
  'location': {'x': -8758500.574766971, 'y': 4276897.574964388},
  'score': 94.81,  'attributes': {'Loc_name': 'World',   'Status': 'M',
  ......

This gets me a geocoded address.   Next, I try to query the feature layer:

requests_layer = FeatureLayer(wake_url)

Setup my SQL statemetn for my where clause:

from_date = dt.datetime.strptime('1/1/2016', '%m/%d/%Y')
sql = "SALE_DATE >= DATE '{}'".format(from_date.strftime('%Y-%m-%d'))
sql

  "SALE_DATE >= DATE '2016-01-01'"

Then perform my queries:

#how many total records are there?
requests_layer.query(return_count_only=True)

  393540

Now, how many parcels have sold in last 4 years:

requests_layer.query(where=sql, return_count_only=True)

  95117

All good so far as I now have filtered this down to just the parcels that sold recently across the whole county.  But I then get lost when I try to further reduce this to just the parcels within a radius of an address...     I have tried buffers, geometry_filter, etc.  Can't seem to find the right combo.  Example:

requests_layer.query(where=sql, return_count_only=True, geometry_filter=address_geo, distance=1, units='miles')

  95117

This is a small subset of the county and the count should be much less.   What am I missing?

I've also played around with the endpoint REST API.  I"m not sure I understand what the geometry_filter is doing for me.


texts: How can i display markers from a database?

texts:
I'm an absolute GIS beginner and here is what i am trying to achieve.
I am developing a web application that allows users to :

type an address, and be displayed the corresponding slipping map
select a point on this map, and create a marker with custom info
this marker and associated info will be stored in a database
every user requesting an address will be presented a slipping map,
showing all the markers other users have defined.

So far, i have learnt how to get lat/lon from an address, using nominatim, and display a map via openlayers. I think i know how to store markers in a database, but i don't have a clue how to display an openlayers slipping map, with my markers (several thousands) applied on it ?
Is it possible without installing my own mapserver/geoserver ?
Thanks for your hints.

Christine


texts:
I need to run a Viewshed analysis (ESRI 10.2 spatial analyst tool) to answer this: what areas will be able to see a proposed 30m high structure.

I have the choice of a DSM or a DTM to use as the input elevation model. I am thinking that the DSM will be better as it will be the tops of trees/buildings etc.) which in the real world will (mostly) obstruct  the views of the structure. Is that the best one to use?

If so what is my Offset B ( the height of the observer). Usually I would put this as 1.5m - does this still make sense to do this. Or will this be like putting a person standing tops of the trees?


texts: Knowing date of BING Satellite images in QGIS QuickMapServices plugin

texts:
I really like using QMS's BING satellite basemaps because of their very high resolution. They are the best I've been able to find. But I need to know when the basemap layer images were taken.
Is there any way to find the date for the Bing Satelite (QMS plugin)images displayed?


texts: Finding boundary between two different groups of points

texts:
I have a dataset with premises and their location in lat/lon points. Not all the premises are the same. Some premises are offices and other premises are houses. I am trying to find the best way to create polygons that cover as many premises that are houses while covering the least amount of offices. For simplicity, in this question we can consider 1 house to be worth the same as 1 office.
I currently have a good working solution by creating a grid of polygons on top of all the premises and calculating how many houses and how many offices are in each cell. Then just grouping the cells that cover more houses than offices into bigger shapes.
However, I am wondering if there is any algorithm that can create a boundary line between the two types of premises? There is no need to divide the premises into 2 different areas, this division can be in many different smaller areas, positive areas within a bigger negative zone etc..
There has to be a better, smarter way to do this. It feels like my current approach is just brute force. Any ideas?
EDIT
I am currently running my approach with Python using GeoPandas and Shapely files


texts: Matplotlib Basemap how to set EPSG to be consistent with new version of pyrpoj

texts:
I am getting a Matplotlib Basemap error related to pyproj and setting the EPSG. (I know Basemap is being depreciated, but this is what I am working in for now.) As a result, I can't plot the arcgisimage basemap.
I'm working in a new env, so I think it is related to a pyproj update (everything worked perfectly in the old environment). I think I need to change how the EPSG is set in m = Basemap() in my code below, but I can't find the correct syntax in my searches.
Error:
C:\myenvpath\lib\site-packages\pyproj\crs\crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
C:\myenvpath\lib\site-packages\pyproj\crs\crs.py:294: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  projstring = _prepare_from_string(" ".join((projstring, projkwargs)))

Code:
south, north = 34.2, 34.5
west, east = -118.3, -118.0
center = [(east + west) / 2, (north + south) / 2]
m = Basemap(llcrnrlon = west, llcrnrlat = south, urcrnrlon = east, urcrnrlat = north,
                resolution = 'c', epsg = 4326, lon_0 = center[0], lat_0 = center[1],
                suppress_ticks = False)
img = m.arcgisimage(service = "World_Shaded_Relief", xpixels = 2000)
img.set_alpha(0.5)

Using Jupyter Notebook & running the following:
Windows 10
conda 4.8.2
Python 3.8.3
basemap 1.3.0 py38hcdd3ad8_2 conda-forge
matplotlib 3.2.1 0 conda-forge
matplotlib-base 3.2.1 py38h1626042_0 conda-forge
pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge


texts: ST_Transform result value of PostGIS is different

texts:
PostGIS is installed on each of the two servers.
I run the following query on one server to get the following results:

query

select st_astext(st_transform(
      ST_GeomFromWKB('01010000203B1400001CB72B9CCEC62D41C14E439A5EBF3D41'::geometry, 5179),3857))

result1

POINT(14162637.4168507 4515322.6138588)

and then I ran the same query on the other server and got the following results

result2

POINT(14162637.4168507 4489271.07636558)

The two servers' spatial_ref_sys tables had the same definitions for 5179 and 3857, respectively.
I do not know why two result values have the same x value but different y values

result1 server version is ...

POSTGIS="2.1.2 r12389" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.6" LIBJSON="UNKNOWN" TOPOLOGY RASTER

result2 server version is ...

POSTGIS="2.4.5dev r16619" PGSQL="96" GEOS="3.6.2-CAPI-1.10.2 4d2925d6" PROJ="Rel. 5.0.0, March 1st, 2018" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.9.1" LIBJSON="0.13.99" RASTER

is there a bug in a particular version?

epsg:5179 definition

5179    EPSG    5179    PROJCS["Korea 2000 / Unified CS",GEOGCS["Korea 2000",DATUM["Geocentric_datum_of_Korea",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6737"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4737"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",38],PARAMETER["central_meridian",127.5],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1000000],PARAMETER["false_northing",2000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","5179"]] +proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs


texts: What does it mean by NULL geometry of a shapefile in ArcGIS?

texts: I was trying to work on a line shapefile and while trying to call the firstpoint of the shapefile, it was throwing me an error saying the shapefile has nothing called firstpoint. Then I repaired the geometry using the ArcGIS tool and it got fixed. But, I want to know, what did that mean by repairing the geometry of the shapefile.


texts:
I just upgraded to ArcGIS server 10.3, and I still can't print with secure map service.

I embedded username and password in my web Application (visual studio .NET 2013 C#).

I saw this issue in the previous versions in the error list. But I thought they fixed this issue in 10.3 .
Does anybody else still have this problem?


texts: Change Detection ArcGIS with multiple dates

texts:
I'll explain my aim first: I am trying to create a single raster to represent vegetation resilience, I am doing this via change detection.

So far I have created 5 Change detection images, where 2000 is a reference image. I have done nothing to these except to display them with the classify tool with 1 st dev.

My problem is what is the best way to combine these images to demonstrate resilience. I'm assuming it's just a raster calculation.


texts: R and osmar : How to merge two adjacent osm files with overlapping nodes and ways?

texts:
Due to a problem loading large osm files in osmar, I have decided to split my osm file in small parts, load them individually and the merge them.

I currently face two problems :

The code splitting the osm file tends to "lose" data around the bounding box. Therefore, I can't merge the objects afterwards because I have lost information in the process. This is done using osmconvert : osmconvert filename.osm -b=4,44,5,45 -o=filename_out.osm for your information.
To solve this, I now take a margin around each bounding box to make sure I don't lose nodes or ways. Consequently, I have to merge objects with overlapping nodes / ways in some cases. This is done with c() applied to objects of class osmar. And the result is not satisfying (see images below, with a map split in 3*3 tiles and then combined again) :  and

The code used to display these maps is spplot(as_sp(osmar_object, "lines"), zcol=4) for your information.

How am I supposed to merge these files ? Should I find a way to merge osm files with overlapping nodes and ways or should I learn a better way to split the file than osmconvert ?

I'm opened to any approach, including switching to sp objects if needed.


texts:
Trying to find the maximum of a several Shapefiles I execute a cascade of unions. Out of all the tools, so far GRASS v.overlay seems to be the one working the best. However I am running into trouble with a dataset and it seems to be a topology issue.

Consider this example:
I have two Shapefiles/layers.

Input 1:

Input 2:

The resulting union:

This first union works as expected, even if the resulting polyong is very small. Long diagonal of the yellow polygon is about 1mm. The field I want to maximise is an Integer called "int" and the resulting table is the upper one:

However if I execute the next Union with another polygon which covers all the area displayed in the images completely with a single variable (so no more slivers, yellow polygon is identical) the table suddenly misses values of the Unions already made. E.g. a_a_int in the second table should be a_int of the first (0).
All layers are in the same coordinate system. This happens only for a few polygons in the whole process and they seem to always be very small. How can I avoid this/work around?

QGIS 2.18.16
GRASS 7.4.0


texts: Showing cumulative load by road segment using QGIS?

texts:
The picture shows the forest roads where wood transport is carried out by trucks. From polygons, logs (wood) are collected near roads and then transported.

How do I show or sum cumulatively  collected the total load by parts of roads using QGIS?

This is necessary so that I can see which section of the road will  be damaged soon. Therefore, it could suggest which part of the road will be reconstructed. Units are in tonnes (t).

In any part of road I can automatically store loads, but I do not know how to sum loads like on picture shown below and put the new value in road segment, automatically.


texts: PostGIS - not a valid Win32 application error

texts:
I'm trying to install PostGIS but with no luck. I've installed PostgreSQL 11 x64 and PostGIS Bundle 2.5.2 x64 and when trying to create an extension in pgAdmin I get:

ERROR: could not load library "C:/Program Files/PostgreSQL/11/lib/rtpostgis-2.5.dll": %1 is not a valid Win32 application. SQL state: XX000

I have Window 10 Pro x64. I have no idea where does this 32 error comes from. The installers come in x64 version so this should be available and I remember having this working correctly on my previous laptop.

I've also tried with PostgreSQL 10 x64 and PostGIS Bundle 2.4.4 x64 with the exact same outcome.

Have anyone experienced this or have an idea on what else to check? I'm stuck.


texts: Using arcpy.SearchCursor() for multivalue field?

texts:
I need to make a multivalue field to be selected as:

Input_Features = arcpy.GetParameterAsText (0)

Num_Pop = arcpy.GetParameterAsText (1) #(It contains population data for the periods 2020, 2025, 2030 and 2035 to about 1800 polygons#

PerPopFic - It contains percentage of data to calculate the new population and overwrite the respective fields

And be calculated for each field chosen Num_Pop a multiplication with PerPopFic field that is in the same shapefile Input_Features.

I've tried some routines, the last I stopped was this:

fc = Input_Features
rows = arcpy.SearchCursor(fc)
arcpy.ListFields fields = (fc, Num_Pop)

for row in rows:
  for field in fields:
     Expres3 = "round ([{0}] * [PerPopFic], 0)". Format (field.name)
     arcpy.CalculateField_management (Saida_Ident_1, field.name, Expres3, "VB", "")

This Num_Pop variable contains on average about 4 fields.

I want to overwrite existing fields. I have 4 fields with population ( 2020, 2025, 2030 and 2035) and another (PerPopFic) with percentage to calculate the new population. All those fields are in the same shapefile.

But I need to do this calculation to another 5 shapefiles who the population field haven't the same name then this. That's why I need to choose these fields.


texts:
Is there a way to subselect features in QGIS as can be done in ArcMap?
For example, I have 5 features selected, and I want to look at each one of them individually (e.g.: zoom to), without unselecting the other 4.
In QGIS when I have a selection, I lose the rest of the selected features if you click on any one of the selected features. I want to be able to subselect (i.e.: create a selection within a selection).


texts: Finding the correct Web Server Base Folder URL - ArcGIS Viewer for Flex

texts:
When getting started with ArcGIS Viewer for Flex I enter the following details;

Web Server Base Folder:
C:\inetpub\wwwroot\flexviewers

Web Server Base Folder URL:
my machine name followed by the \flexviewers (can't provide http address here)

However, I receive the following error message;

Please check your Web Server Base Folder URL.

I can't understand why the URL I provided isn't correct.

I am running on Windows 7 and have turned on II 6 Management Capability and Management Console. Is it possible that I need to share the folder on the server and if so how would I go about doing this?


texts: Convert heights from DHDN Gauß Krüger Zone 4 to UTM 32 ETRS 89 in pyproj

texts:
I am using pyproj.

X/Y conversion is pretty straightforward through

transformer = Transformer.from_crs(31468, 25832,always_xy=True)

However, this does not transform the heights.
From
Vertical Datum transformation using Pyproj
I got the impression that I might need to convert the heights with different CRS since the ones I am using are strict lateral.

I tried initializing the transformer with the corresponding codes for the ellipsoid systems

Transformer.from_crs("EPSG:6314","EPSG:6258")

But that throws

  Invalid projection: EPSG:6314: (Internal Proj Error: proj_create: crs not found)

I need help here.


texts: Popup positioning issue in OpenLayers 2.10

texts:
We were using OpenLayers 2.7 before and everything was working fine.
Then we decided to upgrade our system to use OpenLayers 2.10.

We have popup positioning issue with this new version. The popup is getting displayed away from the actual marker. It is working fine if I roll back to use OpenLayers 2.7.

I can not insert a snap here but you can have a look at the examples below. Ref. http://dev.openlayers.org/releases/OpenLayers-2.10/examples/popupMatrix.html

I guess there is not an integration issue in my system/code as it's working fine with OpenLayers 2.7. Ref. http://dev.openlayers.org/releases/OpenLayers-2.7/examples/popupMatrix.html

I checked it with all browsers, the behavior is same for all.


texts:
If you select in QGIS:  View, Decorations, North arrow, then a window pops up with some parameters regarding the north arrow. One of these parameters is Angle. If you select automatic, then the angle is shown in degrees. Obviously this value is the angle in degrees (in my example 170) between the north position that is currently on the screen and the true north position.

Is it possible to retrieve or to calculate this value (like the north arrow function does)?

I work with WGS84 /EPSG:3031. I have an excel file with ship position data including heading information that I can load into QGIS. It works, every position is marked with a point. But instead of a point I would like to see a little arrow which shows the ship's heading. But this works only when I know this degree value. I could then add the ship's heading information and the arrow would be displayed correctly (with the help of Layer properties, single symbol up-arrow, rotation parameter). In my example the red arrows (the ship's course) are labeled with the timestamp and the heading information from the file.


texts:
When I georeference a raster image in QGIS, it's always a bit displaced as in the attached image (the cyan polygon is the template, the red border is the georeferenced image), although I try to set the georeference points as exact as possible.

Is there a trick to set the georeference points more exact or maybe move the raster layer a bit? Maybe there is a different software, you can georeference by moving and then import it to QGIS?


texts: Alignment/projection issue when opening a shapefile in ArcMap that was sourced from XML file

texts:
I seek to build a UK bus route in ArcMap and am currently  attempting to utilise an .xml file sourced from gov.uk.
I have had some difficulty in converting the .xml file into a shapefile.
So far, I have tried converting the .xml file into a CSV file, adding the X Y coordinates (lat, long) to ArcMap to create a point shapefile. The new shapefile layer was then exported. I added the exported shapefile back into the map and defined its Geographic and Projected coordinate systems and then used the 'Project' tool to attempt to align it with my basemap. These are:

Projected Coordinate System:   British_National_Grid
Projection: Transverse_Mercator
False_Easting:  400000.00000000
False_Northing: -100000.00000000
Central_Meridian:   -2.00000000
Scale_Factor:   0.99960127
Latitude_Of_Origin: 49.00000000
Linear Unit:    Meter

Geographic Coordinate System: GCS_OSGB_1936
Datum:  D_OSGB_1936
Prime Meridian:     Greenwich
Angular Unit:   Degree

These tools failed to align the shapefile and, despite further attempts to ensure the shapefile shares the desired projected and geographic coordinate systems (manual input via ArcCatalog -> Properties -> XY Coordinate System -> New), this shapefile still fails to align with the map.
Can anyone provide any advice on this process? I am new to using .xml files, coordinate systems and GIS in general.
The .XML file can be downloaded from:
https://data.bus-data.dft.gov.uk/category/dataset/152/
I am focusing on the number '46' bus route, which data is included in the JE_46_NMS_PB_46_20200330 .xml file.


texts: Is WFS 2.0 fully backwards compatible with prior versions?

texts: My team has inherited some code that claims to have been written against WFS 2.0 but we need to work with servers that support earlier versions.  Does anyone know whether WFS 2.0 is backward compatible with prior versions?


texts: Control number of Rasters plotted per row in R levelplot

texts:
I want to plot 10 rasters using levelplot but I want to arrange the plots in such a way that there's 1 plot in first row (in the center) and then 3 plots in next three rows. The idea is to show the "control" raster at the top and the month rasters below it (i.e. starting from the next row) so as to make comparisons easier.
#Desired Output looks something like this (numbers being individual plots)
     1
2    3    4
5    6    7
8    9    10

I am using the below code but it plots them sequentially, one after the other (obviously).
library(raster)
library(rasterVis)
library(lattice)
library(grid)
library(RColorBrewer)

#sample data
s <- stack(replicate(10, raster(matrix(runif(100), 10))))
col1 <- colorRampPalette(c("#00008B", "#1874CD", "#00BFFF", "#76EE00", "#006400", "#FFB90F", "#EE7600", "#8B2323"))
mk.attr = c('Control', 'January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September')

levelplot(s, xlab = NULL, ylab = NULL, col.regions = col1, names.attr = mk.attr,
          colorkey = list(space = "right", labels = list(cex=1, font=2, col="brown"),
                          height=1, width=1.4,
                          row=4, column=3, hjust = 2),
          main = list(label = "Main Title", cex = 1.2), margin = FALSE)

#add axis title
grid::grid.text(label = "Trends (in mm/year)", x = unit(0.90, "npc"), y = unit(0.50, "npc"), rot = 90,
                gp = gpar(fontface = "bold", col = "gray33", fontsize = 11))


texts: How to compute geodesic area in GeoDjango?

texts:
I have some polygons which can be anywhere in the world. I would like to compute their areas.

The GeoDjango docs for GEOSGeometry.area don't specify whether it is geometric or geodesic area, but given that GEOSGeometry.distance() says that "GEOS does not perform a spherical calculation even if the SRID specifies a geographic coordinate system", I suspect that it is a simple geometric area.

I'm willing to write a little code, but the shorter the better, of course.

We do have a working geodesic distance function already.


texts: Create a generic shape and draw it as a polygon with center (or base) and rotation?

texts:
In Java (jts), is there a way to draw a generic shape, eg. bounded in meters, and apply this to a coordinate, with rotation around a center or base?

Something like:

double lat = 65.100
double lon = 7.123
Box box = new Box(10, 100, METERS);
Rectangle rect =
    GeometricShapeFactory.createRectangle(box, lat, lon, Math.toRadians(45));

where lat and lon are coordinates in WGS84


texts: Adding Images/Attachments to Survey123 in Inbox

texts:
Project Background: I have a spreadsheet of existing assets that I will be inspecting in an upcoming site visit. I uploaded the Existing Assets CSV file to AGOL to create a feature layer. A PDF of the feature layer screen is here (Lake Pleasant Feature Layer) and the JSON is here (Lake Pleasant JSON). I then created a new survey in Survey123 Connect using the Feature Service radio button, and successfully created a survey in which I'm able to use the Inbox to view existing assets and modify any field I'd like.

Question: I want to be able to take photos of these existing assets during the site visit. I've read online that this is a shortcoming of Survey123 if modifying existing records, and that a repeat field needs to be added to create a related table. I attempted to create a repeat field for the images to be taken in the field; however, I get the error message "The custom feature service submission_url is not compatible with this survey (Table PHOTOS not found)" (also in the screenshot below).

How can I fix this problem?

XLSForm

Existing Assets

Lake Pleasant Feature Layer

Lake Pleasant JSON


texts:
I am trying to follow PostGIS ST_Buffer Radius Help to use a 50 meter distance but postgis is returning an error.

I am trying the following SQL,

SELECT id FROM table
WHERE ST_DWithin(ST_GeogFromText('SRID=4326;POINT(-3.165356 55.926665)'),
      geography(geometry), 50);

The hints indicated at ST_GeogFromText are,

function st_geogfromtext(unknown) does not exist at character 92
No function matches the given name and argument types. You might need to add explicit type casts.

Any help/alternatives would be appreciated.


texts: Set Negative Values to Zero for raster in arcpy

texts:
I have a DEM layer containing floating point data.
I am using the following code to convert the negative values to 0 in arcpy.

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:\mmrahm

out = Con("cnph" < 0.0, 0.0, "cnph")

The code runs fine and a raster layer named "out"  is created. However, the negative values are still there, i.e., the output layers remain identical to the input layer.

Can anyone help me to resolve this issue?


texts:
I'm trying to convert between WSG84 and ED50 (non-UTM) using C# and DotSpatial.Projections:

using DotSpatial.Projections;
using System;
using System.Linq;

namespace ConvertWGS84ToEd50
{
    class Program
    {
        static void Main(string[] args)
        {
            double[] xy = { 44.59740000, 10.69170000 };
            double[] z = { 98 };

            Console.WriteLine($"Input  -> Lat:{xy[0]} Lon:{xy[1]}");

            var wgs84 = KnownCoordinateSystems.Geographic.World.WGS1984;
            var ed50 = KnownCoordinateSystems.Geographic.Europe.EuropeanDatum1950;

            Reproject.ReprojectPoints(xy, z, wgs84, ed50, 0, 1);

            Console.WriteLine($"Output -> Lat:{xy[0]} Lon:{xy[1]}");
            Console.WriteLine("Press any key to continue..");
            Console.ReadKey(true);
        }
    }
}

But it gives me the same output as the input. I'm expecting a small variation in the coordinates like (44.598342, 10.6927). Can someone please suggest me what I'm doing wrong?


texts: Adding multiple divIcons to each point in a GeoJSON - Leaflet

texts:
I am trying to add divIcons (SVGs) as markers to a point layer from a GeoJSON in a webapp using Leaflet.  I need to have different sized markers at different zoom levels.

I have asked related questions here and here in my attempts to simply resize the svg on a map zoom event, but this has proved harder than I thought (and harder than I think it should be). I would like to now try to add multiple markers to each point, and toggle their opacity at different zoom levels, but I don't see a clear way to do this.

Here is the code that constructs the map layer from a GeoJSON, and adds the icon:

geoJSONLayer = L.geoJson(newFeatures, {
    pointToLayer: function(feature, latLng) {
        return new L.Marker(latLng, {
                icon: aDivIcon,
                clickable: false,
                opacity: 1
            }).bindPopup(popupContent, popupOptions),
            new L.Marker(latLng, {
                icon: aDivIcon2,
                clickable: false,
                opacity: 0
            }).bindPopup(popupContent, popupOptions); /// adding a 2nd marker in this way doesn't work
    }
}).addTo(map);

var aDivIcon = L.divIcon({
    "iconSize": [12, 12],
    "iconAnchor": Anchorvar,
    "popupAnchor": Popupvar,
    html: '<svg class="anIcon"><use xlink:href="#marker-anIcon"/></svg>'
});

I tried to swap the css class this way on a zoom event, but this leaves the old and new marker visible:

$(".icon14").attr("class", "icon14 icon15");


texts: How to find out which areas are affected by flooding due to heavy rain?

texts:
I'm working on a little project to visualize affected streets of flooding due to heavy rain. I already got a DEM for my study area and filled the depressions with r.fill.dir.
Since I never worked on hydrological analysis with QGIS, I don't exactly know which workflow/tools I should use to determine the affected area of the flooding due to heavy rain. I know that in reality the affected area is a slope which channels the water downwards. I already could modelling this with r.watershed. But is there a way or a tool to calculate the area of the flooding in relation of the amount of rainfall?


texts: How can I delete renderer from OpenLayers 3?

texts:
If I'm creating new renderer using the following code (taken from http://openlayers.org/en/master/examples/geolocation-orientation.html)  ,

map.beforeRender(function (map, frameState) {
        if (frameState !== null) {

            // use sampling period to get a smooth transition
            var m = frameState.time - deltaMean * 1.5;
            m = Math.max(m, previousM);
            previousM = m;
            // interpolate position along positions LineString
            var c = positions.getCoordinateAtM(m, true);
            var view = frameState.viewState;
            if (c) {
                view.center = getCenterWithHeading(c, -c[2], view.resolution);
                view.rotation = -c[2];
                marker.setPosition(c);
            }
        }
        return true; // Force animation to continue
    });

How can I disable / delete it when I want to stop this rendering ?

****** Solution ********

map.beforeRender(function (map, frameState) {
//console.log (disable_render);
if (disable_render === true) return true;

if (frameState !== null) {
    // use sampling period to get a smooth transition
    var m = frameState.time - deltaMean * 1.5;
    m = Math.max(m, previousM);
    previousM = m;
    // interpolate position along positions LineString
    var c = positions.getCoordinateAtM(m, true);
    var view = frameState.viewState;
    if (c) {
        view.center = getCenterWithHeading(c, -c[2], view.resolution);
        view.rotation = -c[2];
        marker.setPosition(c);
    }
}
return true; // Force animation to continue

});


texts: Computational most efficient way to convert Cartesian to Geodetic Coordinates

texts:
As far as I understand the literature there are several ways to convert a set of geoecentric cartesian coordinates to geodetic coordinates

Bowring, 1976
Borkowski, 1989
Lin & Wang, 1995
Toms, 1996
...

What would be the most efficient way to convert a cartesian coordinates to geodetic coordinates - most efficient meaning fastest (and most direct) when being implemented on a computer?


texts: Automatic correction of many topological errors in shapefile?

texts:
I have a database with topological error. In particular I have an overlay of polygons. I would like to correct these errors.

I thought to use an automatic procedure since I have a lot of errors. I thought to use the v.clean command in grass (I use QGIS).

The command works well but it creates many additional geometries.

Is there a method to avoid this problem?


texts: How to expand a floating point raster by one pixel

texts:
In ArcGIS 10.2.2, I have a floating point raster of oceanographic data that boarders a coastline. Due to the jagged edge of the coast, the raster edge only sometimes overlaps the coastline and leaves gaps of NoData in other places - land shown in grey, NoData cells shown in white.

I want to expand the raster all the way to the land polygon so I can infer oceanographic values for points on the coastline, and it is only one cell short. Since it is a floating point raster and not a categorical raster, the Expand zone function isn't appropriate (there are no zones). Is there an equivalent of Expand for a floating point raster?


texts: How can I get domain values associated with a field using ArcObjects?

texts:
I am trying to access the domain of a selected feature.  I have seen examples for VB.NET and C#.NET, but I cannot access the domain of a specified field.

Dim findName As String = Nothing
Dim findValue As Integer = 21 'the domain value for AGRICULTURE
Dim fieldIndex As Integer = m_Feature.Fields.Field("my_field") 'this has a coded domain

'--Check for my_field
If fieldIndex <> -1 Then
    Dim domain As IDomain = m_Feature.Field(fieldIndex).Domain
    If domain.Type = esriDomainType.esriDTCodedValue Then
        Dim cvDomain As ICodedValueDomain2 = TryCast(domain, ICodedValueDomain2)
        For i = 0 To cvDomain.CodeCount - 1
            If cvDomain.Value(i).Equals(findValue) Then
                findName = cvDomain.Name(i)
                Exit For
            End If
        Next i
    End If

When debugging it hangs up on the line with domain.Type, and it shows that the value of Domain is Nothing.  The documentation shows that IField.Domain should return the default domain as an IDomain object.  I know the field has a domain, but I don't know if it is default and am unable to change it.  Does anyone have any advice on how I should access the domain of a field  (of a selected feature in my case)?


texts: How to remove custom auxiliary field in QGIS 3.2?

texts:
In QGIS 3.2 I have created custom auxiliary field (manually). Then I wanted to remove it but it didn't work with this field, worked only with qgis added ones.

Is it a bug or there is another way?


texts: Identify same values and assign code to them using ArcGIS field calculator?

texts:
There is the attribute table the FID_lote_teste fields and selecao both integer.

The selecao field receive zero (0) if the FID_lote_teste field number does not exist in memory.

For example:

First row: selecao gets 0 (zero) because the number 2 from FID_lote_teste field does not exist in memory;

Second row: selecao gets 0 (zero) because the number 3 from FID_lote_teste field does not exist in memory;

Third row: selection receive 1 (one) because the number 2 from FID_lote_teste field already exists in memory;

(...)

6th row: selecao gets 0 (zero) because the number 7 from FID_lote_teste field does not exist in memory;

7th row: selecao gets 0 (zero) because the number 8 from FID_lote_teste field does not exist in memory;

8th row: selecao receive 1 (one) because the number 7 from FID_lote_teste field already exists in memory;

(...)

How do I fill in the selecao field with the field calculator following this algorithm?


texts:
I installed:

gdal-110-1600-x64-core.msi
GDAL-1.10.0.win-amd64-py2.7.msi
gdal-110-1600-x64-filegdb.msi

from www.gisinternals.com

I can open my gdb on the command line:

C:\data>ogrinfo --version
GDAL 1.10.0, released 2013/04/24

C:\data>ogrinfo data.gdb
INFO: Open of `data.gdb'
  using driver `FileGDB' successful.

In "C:\Program Files\GDAL\gdalplugins" I have a file named ogr_FileGDB.dll

However, when I execute the following code:

os.environ['PATH'] = "C:\\Program Files\\GDAL" + ';' + os.environ['PATH']
os.environ['GDAL_DRIVER_PATH'] = "C:\\Program Files\\GDAL\\gdalplugins"
os.environ['GDAL_DATA'] = "C:\\Program Files\\GDAL\\gdal-data"

print os.path.exists(os.environ['PATH'].split(';')[0])
print os.path.exists(os.environ['GDAL_DRIVER_PATH'])
print os.path.exists(os.environ['GDAL_DATA'])

driver = ogr.GetDriverByName('FileGDB')
print driver

I get the following output:

True
True
True
None

It looks like the FileGDB driver does not work with Python.
Did I miss something to install?

Is it possible to print all available drivers of GDAL in Python?

p.s. A similar questions was asked already: "How to Open a File Geodatabase with ogr? [duplicate]". It was marked as duplicate to "How to access feature classes in file geodatabases with Python and GDAL?". However, I think there ogr.GetDriverByName('FileGDB') clearly works.


texts: Adding GeoPackage to ArcGIS Desktop using ArcObjects

texts:
How do I add a GeoPackage layer to ArcGIS Desktop using ArcObjects code (preferably c#)?
From the ESRI documentation:
The WorkspaceFactoryProgID identifies the ProgID of the WorkspaceFactory object that is used to open that particular workspace.
Here is a list of some of the possible WorkspaceFactoryProgIDs:
esriDataSourcesGDB.AccessWorkspaceFactory
esriDataSourcesFile.ArcInfoWorkspaceFactory
esriDataSourcesFile.CadWorkspaceFactory
esriDataSourcesGDB.FileGDBWorkspaceFactory
esriDataSourcesOleDB.OLEDBWorkspaceFactory
esriDataSourcesFile.PCCoverageWorkspaceFactory
esriDataSourcesRaster.RasterWorkspaceFactory
esriDataSourcesGDB.SdeWorkspaceFactory
esriDataSourcesFile.ShapefileWorkspaceFactory
esriDataSourcesOleDB.TextFileWorkspaceFactory
esriDataSourcesFile.TinWorkspaceFactory
esriDataSourcesFile.VpfWorkspaceFactory

The use of the word "some" makes me think that there might be more factories available.  You can certainly add GeoPackage layers using the AddData button in ArcMap.


texts: Creating one layer through use of SQL queries for land use categories

texts:
I have a map with about 63 grids and trying to determine the land cover user category of the greatest area size in each grid.   The land cover categories have been converted to polygons and I already have attribute tables with a field for "Area" in each grid.  I need to create a layer that will give me which land use category has the greatest area size.  So I need to sum up the "area" of each category and then determine the max (highest area) and then make a layer for only that category.

i.e.
Grid 1 contains - water = 900m,  development = 1000m, forest = 2000m
Grid 1 - category with highest area size = forest at 2000m - I would need just one layer displaying forest only in that grid.  But then do this quickly for all other grids in my map.

I am a mid-level user in ArcMap with decent SQL query skills, so I feel most comfortable in SQL and I'm certain a making a query with a join statement that joins onto itself can get that done for me.  I'm looking where/how best to run a SQL query that will be a little more complex than the regular "select * from".  Looking at tools like the "Select attribute by" is nice but the SQL query section in there doesn't really allow for inner join statements to be performed.  Or if all this needs to be done a different way how I go about to create the layers I need. I could of course trudge through each grid and do it manually but I am looking to learn how to do such things in a faster method and use in the future.

Will I be venturing into ArcPy and/or ModelBuilder?

I want to sum up "Area" per "Gridcode" then determine which "Gridcode" has the highest sum area.  Then make that into a new layer.  I would not need any other gridcodes if they do not sum up as the max.


texts: Polylines rotation through the model builder

texts:
I am searching to add rotation to polylines through the model builder.
Is anybody familiar with such option?


texts: Minimum data needed to perform nearby radius query with distance in PostGIS?

texts:
I am new to PostGIS and and PostgreSQL as well. I have been experimenting PostGIS for about a week now and I have seen that almost all blog posts or tutorials require having a shapefile to start querying against.

I can see that PostGIS is really powerful but little overwhelming for me as a beginner since I haven't understand it that much yet.

The main reason I am using PostGIS is to find users within radius and distance calculation for each point found within that radius.


texts: Editing WFS from GeoServer using OpenLayers

texts:
After reading a very good tutorial on how to edit WFS with OpenLayers, I've tried replicating it but with my own WFS layer from Geoserver.

I managed to load the WFS and my basemap successfully and managed to get the buttons to show up. The buttons appear correctly like in the working example from that page but, for some reason they don't work at all even though I copied the entire code.

The code for the WFS to load is:

var geojsonFormat = new ol.format.GeoJSON();

var vectorSource = new ol.source.Vector({
  loader: function(extent) {
   var url1 =   'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=BFTchambers:bft&' +
    'outputFormat=text/javascript&format_options=callback:loadFeatures' +
    '&srsname=EPSG:27700&bbox=' + extent.join(',') + ',EPSG:27700';

// use jsonp: false to prevent jQuery from adding the "callback"
// parameter to the URL
    $.ajax({url: url1, dataType: 'jsonp', jsonp: false});
  },
  strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
    maxZoom: 20
  }))
});

The whole thing together looks like this with the buttons having proper hover, but no functionality.

Also even though my WFS appears correctly along the river Thames when I load it in QGIS, in my app it appears somewhere else in the ocean even though I specified EPSG 27700 (though this is just a minor annoyance at the moment).

My main problem is how to make these OpenLayers edit buttons work?

Fiddle with the whole code (apologies if it looks messy, most of it comes from the working example 2 )
https://jsfiddle.net/Luffydude/ex06jr1e/1/


texts: Selecting points in polygon in a range of time (ArcGIS)

texts:
I have:
 - point GPS for people (time in seconds)
 - buffer GPS for a car (time in seconds)

What I want to do is for each buffer GPS coordinate, I want to find the number of GPS people points inside of it. The only thing is this is time sensitive. I want to have the Buffer GPS time to be 'set' as the time(range of 10 seconds or so). So each Buffer GPS time, only the people points within the specified time range are selected, and this iterates throughout the entire layer.

Unsure if this is even possible with ArcGIS or ArcPy, but if there are alternative options that could be used to accomplish it would be a great help!


texts: Animated lines and line styles in Javascript

texts:
I'm using Leaflet library and Mapbox. I need to create animated line from one city to another, slow motion "travelling" from one point to another point. Is that possible in javascript and how? I can add any tools if it's necessary.

Also, I would need some zig zag line, not just straight line, do you have some advices regards to that? Is there some library or something similar where we can choose line styles maybe?

Thank you.


texts: How to update a field with a count of rows from another table on PostGIS?

texts:
So for a full 2 years I used to simply just write open pg version 9 something like

update tableA a
set columnA = (select sum(column_that_I_need_to_count) from tableB b where st_contains(a.geom, b.geom))

But for some reason updating like this does not work on my new database version 10. It just runs the update indefinitely and never finishes

Now I have to write it as

 update tableA a
 set columnA = sum(column_that_I_need_to_count)
 from tableB b where st_contains(a.geom, b.geom)

Problem is that it says it doesn't allow count on an update

Create table statements are:

create tableA (
id integer,
name text,
geom geometry(polygon, 27700)
)

create tableB(
fid bigint,
column_that_I_need_to_count integer,
name text,
geom geometry(point, 27700)
)

Is there a reason why different PostgreSQL have different subquery rules for updating?

Also how to do the actual update?


texts: Combine iterations a loop as a single raster in R

texts:
I am looking to combine all iterations of a loop into a single raster using R. I am performing shadow analysis where I perform the calculation for specific hours of a day and would like each calculation to be combined together into a single raster. For example when I have five hours to calculate I would like my raster to be: outRast = (i1 + i2 + i3 + i4 + i5).

I am having trouble adding the iterations of the loop to a single raster where they will be combined, rather than it being over written.

Below is my code with the output iterations being saved out as a tif.

#Spring parametres
yy <- 2019
mm <- 03
dd <- 20
hh <- seq(10, 14)

#Shadow Analysis for spring season
for (i in hh) {
  jd <- JDymd(yy, mm, dd, i)
  sv <- sunvector(jd, rastCentre(ffLidar)[2], rastCentre(ffLidar)[1], -4)
  springShadows <- doshade(ffLidar, sv, 0.5)
  writeRaster(springShadows, paste("springShadows", i, ".tif"))
}


texts:
How can I normalize pixel values from 0 to 1 or 0 to 100?
Is there any function or command related to it? I'm using JavaScript.


texts:
I have 207 geodatabases, each with 3 feature classes, and multiple shapefiles (point, line, and polygon) within each feature class. The data is from The National Map (USGS). In the screenshot, you can see three of the geodatabases - Cornwall, Danby, and Dover Plains - with their feature classes. My goal is to merge (or append) these 207 geodatabases into one geodatabase that maintains the same feature classes (most important: foundation, gaz vector, and hydrography). I am not well-versed in Model Builder or Arcpy. A cumbersome method would be to append each feature class individually within each geodatabase with the master geodatabase as the target.

How can I make this process more efficient?