texts: Break a bounding box into smaller boxes (coordinates)

texts:
My objective is to create a collection of bounding box coordinates that I can supply to a search function. The idea is that I have the search function search smaller geographic areas, instead of one large one.

To achieve this, I would like to develop some function that takes the following inputs:

Overall bounding box (to be divided) (4 lat/lng points)
Area of smaller bounding box

Are there any tools that I could use to do this?

The output should be a collection of sets (each has 4 points) representing bounding boxes that are inside the overall (parent) bounding box.

I am happy to write a tool using any language or alternatively use an intuitive web based tool.

Any ideas on how to do this?

I am in the process of figuring this out using simple functions that calculate the distance between two points, but writing here for sanity check also.


texts:
I used the example of this website : http://openlayers.org/en/v3.8.2/examples/feature-animation.html for make an animation.

I don't have problem for stopping the event with the function clearInterval(interval); if I declare before my interval var interval = window.setInterval(addRandomFeature,1000);

The problem is I want to clear the layer who contain the object (after my function clearInterval(interval);) and I don't understand what's the layer who contain the event and I don't know how can I clear this layer.


texts: What are the options to auto populate metadata for data received from external sources?

texts: We receive quite a many spatial data from various sources and export to our enterprise database (Esri). But no metadata is available in Description table. I am looking for a tool or method to get all the attributes, which can be easily populated from the files e.g. spatial extent, projection etc.


texts: Making a 3D polygon with 3D view in qgis 3.4

texts:
I tried to make a 3D map in QGIS 3.4, using the "3D view" function.

However, as I change the 3D property the polygon have not changed.
I tried this with a basic ESRI Shapefile that I made.
I also change the crs to "Sphere_Cylindrical_Equal_Area" thinking maybe the CRS need to be in meters and not degrees, but it did not help.

What is needed to make the layer display into a 3D display?


texts: Assigning value to point feature class based on count of cells with certain raster values?

texts:
I have a point feature class that I would like to assign a new field based on the neighborhood of a raster layer.

I'd like to count the number of cells within that neighborhood and write that count an attribute field for each point feature.

The points are often close enough to where there will be overlapping neighborhoods.


texts:
I have a grid of lat/lon points all over the state of Texas and I'm trying to find the flood zone and base flood elevation for each of the points. I've found the flood zone fine by doing a spatial join, but the base flood elevation layer is line data (cross sections) and I don't know what to do - looking around it seems like I might have to generate a DEM or something? How do I interpolate and find out what the BFE is supposed to be? I just don't know what the process should be.

I have a list of lat/longs (point geometries) and I need to get the base flood elevation for them. FEMA provides cross sections for base flood elevations (i.e. not a raster DEM, but line geometries) and what I need is to convert the cross sections to something I can run an algorithm with the lat/longs as the input and the base flood elevation as an output (added to the attributes table for the points) - or maybe there is some method to run this as-is using the line geometries?


texts: Add variables to an existing netCDF file in R

texts:
I have a netCDF file with one variable 'gpwcount' and its lat/lot coordinates also available. The dimension of this 'gpwcount' variable is 1050 rows x 1050 columns. This file consists of population count data for the year 2015 from three west African countries that had a major Ebola outbreak. The file is attached here

I am trying to model an epidemic and run a simulation. Is there a way to add the following variables to the existing netCDF file using R (or any other tool).

The variables I would like to add are labelled,

 Exposed (1050 x 1050 grid all filled with 0);
 Infected (1050 x 1050 grid all filled with 0);
 Hospitalised (1050 x 1050 grid all filled with 0);
 Recovered (1050 x 1050 grid all filled with 0);
 Dead (1050 x 1050 grid all filled with 0);
 Prevalence (1050 x 1050 grid all filled with 0);
 Inhabitable (1050 x 1050 grid all filled with either 0 if the gpwcount in that cell is 0 or a 1 if gpwcount in that cell is greater than 0);


texts: Adding a legend to points in the tmap package

texts:
I am trying to create a map with several points (about 20), where each point will be represented by a number, and where a legend will have these numbers and names of my points.
Is it possible to do it with the tmap package (or any other R package)?

This is what I have:

library(sf)
#> Linking to GEOS 3.5.0, GDAL 2.1.2, proj.4 4.9.2
library(tmap)
data(Europe)

df_obj <- data.frame(id=c(1,2), name=c("First", "Second"), x=c(0, 25), y=c(45, 50))
sf_obj = st_as_sf(df_obj, coords = c("x", "y"), crs = 4326)

tm_shape(st_as_sf(Europe)) +
  tm_polygons() +
  tm_text("iso_a3", size="AREA", root=4) +
  tm_shape(sf_obj) +
  tm_symbols() +
  tm_text("id", col="name") +
  tm_format_Europe() +
  tm_style_bw()

This is what I want (see the legend):


texts: Download data from neighbourhood with OpenStreetMap

texts:
I am wondering how I can download all the streets, and their housenumbers from a neighbourhood from OpenStreetMap. The neighbourhood that I am interested in is located in Rotterdam. I want to import the data into a Neo4j graph database.

How can I download all the streets, and their housenumber from OSM? And what is the best format to import it into Neo4j.


texts: For some reason the Main Menu dis-attached itself from the ArcCatalog window and I am not able to dock it or even move it when i click the menu title bar. The functions work but it's just stuck in place on screen and I can't do anything with it.  The title bar itself appears dimmed.


texts: Using CartoDB for citizen comment collection?

texts:
I'm wondering what would considerations I should make for using CartoDB as a tool for collecting geo-located citizen comments.

Currently we're using a javascript / bootstrap application that uses a google map to get a users X,Y location. It requires the manual entry of an address or the movement of the pin on the map to get that X,Y. There are then other form inputs where the user can enter their name, and other comments about the location, their contact info, etc.

The form inputs, as well as the XY, are then fed to a google spreadsheet where they are manually extracted to a CartoDB map on our main site.

Can CartoDB be used to replace this google doc method, and what do I need to think about in terms of using the SQL API?

Is it a good practice to expose the necessary information to use that SQL API?

How could I use the 'get my location' feature on most browsers to automatically get a users location and feed it into the geometry of row in a CartoDB table?


texts: Brackets meaning on EE visParams of Map.addLayer?

texts:
In the context of https://docs.google.com/document/d/1mNIRB90jwLuASO1JYas1kuOXCLbOoy1Z4NlV1qIXM10/edit# section 3, there is one commented line:

// Map.addLayer(trend, {}, 'trend array image');

From Map.addLayer(eeObject, visParams, name, shown, opacity) documentation, I see it refers to:

visParams (FeatureVisualizationParameters|ImageVisualizationParameters, optional):
The visualization parameters. For Images and ImageCollection, see ee.data.getMapId for valid parameters. For Features and FeatureCollections, the only supported key is "color", as a CSS 3.0 color string or a hex string in "RRGGBB" format.

Is that right? May I ask for a more elaborated explanation of it?


texts: An Overpass query for OSM ways; both part of a relation and also a specific tag?

texts:
How do I construct an Overpass query to find data on Openstreetmap which is both part of a relation and has itself a specific attribute?

Specifically: I'm wanting to find sections of cycle route (in a route relation) AND which are surfaced a particular way using Overpass - accessed using the Quick OSM plugin on QGIS.

So in terms of the OSM data I'm looking for a relation with type=route and route=bicycle
http://wiki.openstreetmap.org/wiki/Cycle_routes#Relations

But I want to limit this by surface - lets say surface=paved or surface=asphalt
http://wiki.openstreetmap.org/wiki/Key:surface

Ideally I'd like to do this using the QGIS Quick OSM plugin - but if necessary I'm happy to use the Overpass Turbo website.

I'm well out of my depth on getting the syntax correct for this (although I've used all these systems for simpler queries).


texts: For GeoServer, it is necessary to create a new style for each different layer to load?


texts: How to calculate distance from point to polygon with a time frame?

texts: Using the Near tool, I knew that I can just calculate the distance of each house to the nearest fire. However, I only want to keep fires that happened in the past 7 years for each house. For example, I have fire data from 1950 to 2013, a house was built in 2002. I only want to measure distance to the nearest fire from 1995 to 2002. Any idea on how to do that?


texts:
I am trying to use the expression editor on the legend item in the layout manager. My layer has the attributes of RU_NAME which is the geological unit and LITH_SUMM which is the description. Unfortunately the layer is huge and has colouring and categories set up by our local government department so I cannot just reclassify without loosing all of the colour and symbology styles.

Concat(@symbol_label,': ',"LITH_SUMM")

Above is the code I am using which should give me the combined name and description of the rock units. The preview has the correct output, but when I click OK it results in [%Concat(@symbol_label,': ',"LITH_SUMM")%]

Is there something I am doing incorrectly with the code for this expression?

EDIT: Please see the screenshot below:

Our local government has released the shapefiles and styles for geological layers which are compliant to standards so I would not want to recategorize the layers based upon a new attribute. Doing this would overprint the colours and labels on a data set that has thousands of entries. The issue is that the legend only shows the symbols "Tp" instead of the description "Arkose sandstone with interbedded siltstone". The information is in the attribute table of the layer, I am just unsure how to get it to be displayed in the legend.

EDIT 2: See Below for additional information

I have attached an excerpt of the attribute table and the formatting style. The symbology has been set by the government to a standard colouring based upon RU_NAME. There are thousands of these categories. I have combined the required information in the attribute table asLegendText = Concat("MAP_SYMBOL","RU_NAME","LITH_SUMM") as an example of what I am trying to get to display in the legend text. I have tried to edit the layer style python and this guide:Using attribute data for legend labeling in QGIS? but I am definitely a novice at python and cannot get it to work.


texts: Using sDNA to analyze the network downloaded from OSM, do we need to merge the network first？

texts: In the road network downloaded from OSM, each road is divided into small sections. I would like to know whether I need to merge the sections of the road network of the same road into a complete road before using SDNA to analyze the road network. What tools to use if necessary?


texts: Summarizing values per row using field calculator via ModelBuilder?

texts:
I think my problem is easy to solve, but so far I couldn't find a solution. I have a Feature Dataset created with the Union tool in ModelBuilder and now try to sum up all values of one specific field (VSG_EINZEL->up to 150 fields) into a new field (VSG_ALL) automatically in ModelBuilder.  My problem is, that the number of VSG_EINZEL-fields is rather high and varies so I can’t use a static field calculaion like:

VSG_ALL = VSG_EINZEL + VSG_EINZEL_1 + VSG_EINZEL_12

and so on.

I am very new in using Python.


texts: Create New Feature Class within ArcCatalog

texts: I am trying to create a new feature class and when I right click on my project in Catalog, I do not have the option for NEW.  See screenshot.  Is this a permissions issue? Is the software not setup properly?


texts:
I have a fishnet layer where some empty cells are nested within. I'd like to only show the boundary of the filled cells (color set no fill, but no inner cell boundaries of adjacent cells). How can I best delete the boundaries between adjacent cells?

I'm using ArcGIS Pro


texts:
Aa while ago I drafted a semi-universal guide on how to compile GRASS GIS on Ubuntu.

Because I had enough emails asking more details I thought to post it here for everyone and if possible to evolve it as we go!


texts: QGIS for mac acting strange after upgrade to 2.16

texts:
I'm running OSX 10.10.5 and just upgraded QGIS from 2.14 to 2.16.

1st. the installed GDAL 2.1 wasn't registering when using gdalinfo --version. It still showed 1.11... So I used MAPBOX GUIDE to get the 2.1 to be added to my PATH. But running SAGA tools from Processing seems to be broken ex. running Fill Sink only results in ERROR.

I tried to re install the 2.14 again (removing qgis.app and ~/.qgis2)

Now 2.14 which worked earlier is not working properly either!

Could anyone give me a clue how to properly install either 2.14 or 2.16 on a mac so both SAGA tools and specially Raster tools will work.


texts: Squared buffered polygon based on polyline

texts:
I'm using ArcGIS Javascript API to allow users to draw polylines and polygons onto a map. Now when a user draws a polyline, this should result in a polygon surrounding the polyline at a fixed width:

(polygon here is handdrawn, but should give an idea of what I would like to accomplish)
I can utilize the GeometryService.buffer method, but that gives me a rounded polygon:

Is there any way I could achieve this? It doesn't seem the BufferParameters object has an option for that. I was also considering using a convex hull might be a way towards a solution, although I'm not seeing yet how I could utilize it to my needs...
Does anyone have a solution?
Update
For now I have achieved a workable yet not ideal solution:

Take the array of coordinates ('rings') from the buffered polygon geometry
Loop through the array of coordinates and calculate the distance between them
If the distance between two vertices is below a certain treshold then I know the coordinates are part of an arc.
In my case an arc of 180 degrees consists of 49 rings, so I remove all rings from the arc except the 25th.
Restore the buffered polygon geometry with the new, reduced set of rings

It doesn't produce the exact shape I'm after but for the time being this is sufficient:


texts: CREATE Script for Enterprise Geodatabase?

texts:
I have an ArcGIS For Desktop license, and we have a PostgreSQL 9.2/PostGIS 2 database server running in our network.  We currently use Query Layers to add the PostgreSQL data to ArcMap sessions.

We would like to make full use of ESRI's geodatabase functionality, without committing to the (we have heard) somewhat unreliable FGDB format.

I see that using the "Create Enterprise Geodatabase" geoprocessing tool one can create a database and have it set-up as a geodatabase, but this requires the input of an authorization key, which it seems is only available with a Server license (true?).

My question:
Can we not just turn this diagram, or better, this schema into a series of CREATE statements for the PostgreSQL database in order to generate the equivalent database schema ourselves, and then utilize it from ArcCatalog as we would a file geodatabase, without any "authorization key"?

My thanks for the clarification.


texts: ArcGIS hyperlink script to assesor maps based on different account numbers

texts:
I am adding a hyperlink within the attribution to open up to our local county assessor maps. I am able to get the hyperlink to work directing the user to the map directory to look up based on a given account number. However, I wish to write a script to use the account number within the attribution to automatically generate the account number's map by-passing the account search page on the website?

So I am trying to hyperlink "ASSESS_LNK" to open the corresponding "ACCT_NUMB" page?

I am using ArcMap desktop 10.8.  So far I am just using a simple 'populate field with text' that points towards the url within field calculator.


texts: Interactively deleting tiles of aerial photographs using ArcGIS Desktop?

texts:
I have a bunch of separated aerial photographs (res. 20 cm, count ~ 50, total size 20GB) stored in single folder. I do not have a reference polygon layer noting the extent and name of each title raster.

My study area is however covered by only part of them. Thus, I wish to delete the tiles (unique rasters) I don't need from my home folder.

Is there a way in ArcGIS 10.1 how to interactively delete the tiles which I don't use? i.e., For shapefiles it is pretty straightforward: in editing mode,
-> select the shapefile -> select a polygons -> delete the unnecessary ones.

I know that I can remove the tiles manually one by one, but how can I make it easier?

Do I need to firstly create a mosaics/geodatabase of all aerial photography and then delete the unused ones?

I seek advice/tips or tutorials about basic managing of tiles of aerial photography.

As an output, I wish to obtain a set of separated aerial photograph - or, ideally - a mosaics of raster images within the extent of the overlapping polygon. I am not sure which one is easier to obtain and display.


texts: Inputting CSV in QGIS graphical modeler and defining its geometry

texts:
Is it possible to have a model in QGIS graphical modeller which accepts a CSV with geometry fields and then allows you to use them for the rest of the model as if it was a vector file?

This is required so I can bypass having to import the CSV into QGIS and defining it there before adding it to a model.


texts: Ensuring Maplex label engine is enabled on new layers?

texts: I have recently run into problems with maplex label engine not being enabled on new layers.  I played around in the Label Manager and that didn't work immediately.  It randomly started working on the new layers 10 minutes later.  I would like to know how to make sure maplex label engine is enabled for any new layers I add to future maps so I can place the labels properly.  It was a Placement Properties issue I was having, only showing the default options instead of the maplex options.


texts: GeoDjango with PostGIS: distance_lte filter is inaccurate

texts:
My view code (relevant part):

location = self.GOOGLE_MAPS.geocode(lookup)
if len(location):
    page = int(request.POST.get('page')) or 0
    start = page * settings.PHOTOS_PER_BATCH
    end = (page + 1) * settings.PHOTOS_PER_BATCH
    geo_range = int(request.POST.get('range',  self.DEFAULT_RANGE))
    geo = location[0].get('geometry', {}).get('location')
    pnt_str = 'SRID=4326;POINT({} {})'.format(geo['lng'], geo['lat'])
    pnt = GEOSGeometry(pnt_str)

    if geo_range > 0:
        # short version
        # photos = Photo.objects.filter(photographer__locations__point__distance_lte=(pnt, D(km=geo_range)))

        #debugable version
        locations = Location.objects.filter(point__distance_lte=(pnt, D(km=geo_range)))\
            .distance(pnt).order_by('distance')
        photographers = Photographer.objects.filter(locations__in=locations)
        photos = Photo.objects.filter(photographer__in=photographers).order_by('?')
    else:
        photos = Photo.objects.all()

    #render the view
    context = {
        'formatted_address': location[0].get('formatted_address'),
        'photos': photos[start:end],
        'geo': geo,
        'total': len(photos)
    }

    return render(request, self.TEMPLATE_NAME, context)

Distance display is made by a custom django filter:

from django import template
from proto.models import Location

register = template.Library()

def distance(origin, destination):
    return int(destination.distance(origin) * Location.SCALE)

register.filter('distance', distance)

Finally, the template code:

{{ point|distance:p.locations.all.0.point }} KM

The problem is that even when the range is set to 1000 M (10 KM), results that seem to have the distance of 14 of 16 KM are still included.
With a range of 20 KM, result of 25 and 28 (and so on) are also included (although they are excluded from the previous result set).

Is it and SRID Issue?

Also, I don't quite understand why should I multiply by 1000 to get the filtering done, and then divide by only 100 to display the distance.


texts: Deleting duplicate points except those with different attributes in QGIS?

texts:
I have a question regarding removal of duplicate points, which I did not find answered here yet. In my point shapefile there are some points at the same coordinates. Some of them have the exact same attributes and some have different attributes.

Is there a way to identify and delete the "real" duplicates (those with the same coordinates AND attributes) leaving out those which have the same coordinates but differ in the other attributes?

I only found tools yet, which only go by the coordinates in identifying duplicates but ignore the attributes.


texts: QGIS capture screenshot of visible map (PyQGIS or other alternatives)

texts:
I would like to take a screenshot of my QGIS area, which I am working on.
It does look pretty much as you can see below:

I would like to grab only the visible area of my map. I am not talking about map canvas, as mentioned here:
Save map canvas as image with Python in QGIS
PyQGIS make screenshot of mapCanvas after setExtent is called
or here
Saving Map Canvas as PNG with transparent background programmatically with QGIS?
or here
Export Map Canvas QGIS
I simply want to have an option to grab the screen from my QGIS working area defined by it visible bounds (after dock windows visible).
I tried some of the codes, but they didn't work:
 mainPath = 'C:\\Users\\mk\\Desktop\\AX\\Pre survey\\Tyton\\'
 filename = 'StrandAAA'
 imageType = "png"

 imageWidth_mm = 400
 imageHeight_mm = 160
 dpi = 300

 mapRenderer = iface.mapCanvas().mapRenderer()

 c = QgsComposition(mapRenderer)
 c.setPlotStyle(QgsComposition.Print)

 c.setPaperSize(400, 160)

 x, y = 0, 0
 w, h = c.paperWidth(), c.paperHeight()
 composerMap = QgsComposerMap(c, x ,y, w, h)
 c.addItem(composerMap)

 # get all items and disable the backgrounds
 itemList = c.items()
 c.removeItem(itemList[2])
 c.removeItem(itemList[2])
 itemList[2].setBackgroundEnabled(False)
 itemList[3].setBackgroundEnabled(False)
 # add them to the composition
 c.addItem(itemList[2])
 c.addItem(itemList[3])

 c.refreshItems()
 c.refreshDataDefinedProperty(QgsComposerObject.AllProperties)

 c.setPrintResolution(dpi)
 dpmm = dpi / 25.4
 width = int(dpmm * c.paperWidth())
 height = int(dpmm * c.paperHeight())

 # create output image and initialize it
 # image = QImage(QSize(width, height), QImage.Format_ARGB32)
 image = QImage(QSize(width, height), QImage.Format_ARGB32_Premultiplied)
 image.setDotsPerMeterX(dpmm * 1000)
 image.setDotsPerMeterY(dpmm * 1000)
 #image.fill(0)

 imagePainter = QPainter(image)

 # This part does not seem to have any effect
 #brush = QBrush(Qt.BrushStyle(0)) # transparent brush
 #
 #imagePainter.setBackgroundMode(0)
 #c.setBackgroundBrush(brush)

 c.renderPage( imagePainter, 0 )
 imagePainter.end()

 imageFilename =  mainPath + filename + '.' + imageType
 image.save(imageFilename, imageType)
 print 'image saved'
 print 'done

but it didn't work, the same as:
  vlayer.select(ID)
 qgis.utils.iface.actionZoomToSelected().trigger()
 vlayer.deselect(ID)
 qgis.utils.iface.mapCanvas().saveAsImage('C:\\Users\\mk\\Desktop\\AX\\Pre survey\\Tyton\\feature.png')

The solution is shown on the YT channel here also wasn't successful
https://www.youtube.com/watch?v=17AXqO5_s7g&ab_channel=OpenSourceOptions
I think I am closer to this solution:
Is it possible to take automatically a "screenshot" of my selected feature in QGIS using Python?
butt instead of selected feature I need the visible map area
The typical solution for taking a screenshot in Windows doesn't solve my request
https://stackoverflow.com/questions/2846947/get-screenshot-on-windows-with-python
although in Pyscreenshot v2.3 I can define the rough bounds of my screen and capture them, how about if I change my QGIS window to another place
https://pypi.org/project/pyscreenshot/#files
Is there an option to capture the screenshot of the visible map in QGIS both by using PyGIS or maybe not programmatically (GRASS or sth?)


texts: How to open memory data with ogr in rasterstats function?

texts:
I would like to reuse vector data created in memory in another function called rasterstats (source code, documentation):

I have a function that creates vector data:

def createPolygons:
    memoryDriver = ogr.GetDriverByName('MEMORY')
    memoryDatasource = memoryDriver.CreateDataSource('memoryData')
    memoryLayer = memoryDatasource.CreateLayer('memoryLayer', geom_type = ogr.wkbPolygon)
    [...]
    return memoryDatasource

I would then like to use the output (polygons) of that function in:

zonal_stats(vector, raster, band=1)

Which results in errors like this:

ValueError("Can't parse %s as a geojson Feature object" % obj)

ogr.Open() in rasterstats seems unable to open my vector data. I have tried returning the memoryLayer and the memoryDatasource in createPolygons and using those as well as their names as input in zonal_stats. The memory data certainly exists, because I can access it after debugging the result of createPolygons with functions such as GetFeatureCount().

The problem is that the data does not have a path, and this seems to be required as input for ogr.Open(). However, memory is listed as a possible driver for ogr.Open(). This side says also: There is no way to open an existing Memory datastore.

I am unsure, whether it is possible to reuse the memory data in another function and if yes, how to realize this.


texts: How to change color ramp interval size in QGIS 2.12?

texts:
I'm working on elevation data and I'm trying to change the interval size to 50 meters in the layer's style. I created a color ramp color interpolation on the elevation layer. The values are between 8-575 meters and the number of the intervals are 12. All the values are integers.

The problem is that the interval size isn't 50 meters, but 51.55
meters and I couldn't find an option to do define it at exactly 50 meters.
Why are the values in decimals instead of integers (50 meters)?

In arcGIS there's an option to change the interval size (also, in the same conditions- number of interval is 12) and the values are integers:


texts: Finding external boundary of bunch of vectors (using R libraries)

texts:
I have vectorized raster as in the image. This is actually 1419 polygons in SpatialPolygons class object.

I'm looking for a method which would allow me to extract the one most external polygon (boundary shape). I've found concave method which works for 95% of cases but not in this particular one. As follows the code:

sfn <- sf::st_as_sf(the_polygon)
sfu = sf::st_union(sfn)
gpts = sf::st_coordinates(sfu)
gmat <- matrix(c(gpts[,1], gpts[,2]), ncol = 2)

I can adjust parameters:

cm = concaveman::concaveman(gmat, 0.5, 1)
chp = SpatialPolygons(list(Polygons(list(Polygon(cm)), ID=n)))

this gives such vector:

but choosing the first parameter as 0.1 takes:

So there's no optimal shape, I've tested values between those 2 and further as well. The second parameter doesn't make much - it just smoothes the lines.

I'd like to ask from your experience what kind of approach I'd take? For example, if I had a function which fills the outside area, then selects "negative" and makes borderline and that border would be good approximate of what I need. Notice it's a vector(s). Have anyone some idea? In terms of fill the vectors area what libraries I'd choose?

  EDIT

Incredible things happen: I've fit this almost ideally using concaveman(gmat, 0.01, 2.2) - just tested with "bruteforce plot". So, in this case, answer been partially found.


texts: How to use offline MbTiles with Telerik's WPF map control?

texts:
Telerik's UI for WPF includes a map control: UI for WPF - RadMap.

There's some discussion of how to configure it to use a local mapping service:
Implement custom map provider, but it seems to be oriented towards directories of tiles: Offline Map Server.

Anyone have any ideas on how to make this work with an MbTiles file?


texts:
New to gvSIG and trying to get SQL filters to work. I can get labels to draw for all features the same way, but when I try to define label classes with filters nothing will draw. I'm trying to label street centerline data based upon a roadway type (Interstate, local, etc)

My condition field is an integer and the filters are:

Update: I've tried joining my numeric values to a table to bring text values over (i.e. my filter now reads 'dm_description='Interstate',etc. Still, the label classes are not being honored.

I'm stuck.


texts: Import the EPSG:7775 projection without PROJ.4 to QGIS?

texts:
I wanted to define the EPSG:7775 projection in QGIS but it has no field for PROJ.4. Everywhere I've seen that you can import a projection with PROJ.4. How do I import this?

I am trying to do this (some basic spatial analysis) for the area specified (Uttar Pradesh) in EPSG:7775.


texts: An old point layer created from PLSS grid does not line up with current PLSS grid layer

texts:
I have two vector layers which I am having trouble lining up. First, I have a point layer that was developed a decade ago to sample information from Public Land Survey System records in Washington and Oregon State. Second layer I have is PLSS grid polygon layer recently downloaded from Bureau of Land Management. The points are not just nodes, but also were placed at various distances along the grid lines.

Spatial reference:

Point layer is projected in NAD83 / BC Albers, EPSG: 3005.
PLSS grid polygon layer is not projected, instead it has GCS North American 1983, EPSG: 4269

I have tried to project the grid layer to EPSG: 3005 and several others. Similarly, I have tried to reproject the point layer. Essentially, sampled points should line up with PLSS grid corners and other point placed along the grid lines, however all of my attempts to project each layer have yielded no desirable results, points are always off. I have tried OTF projection variations in QGIS and project function in ArcGIS Pro. Unfortunately, no documentation is available on development of the point layer, therefore I do not know whether previous PLSS layer used a different projection, but I would expect BLM PLSS layers to be the same to layer available a decade ago.

Please help to identify possible solutions, as the point layer contains a lot of valuable information which would take very long to replace or correct manually.


texts: If/Else in ModelBuilder not showing the correct True/False outcome

texts:
I'm creating a model that queries the number of features in a feature class and runs one of two scripts depending on if the feature class is empty or not.

I'm using Get Count to calculate the number of features with the output being 'Feature Count'

I then have two calculate values 'Policy Detected' and 'No Policy Detected' both set to Boolean data types

Policy detected is:

FeatureCount(%Feature Count%)
def FeatureCount(n):
if n> 0:
    return True
else:
    return False

No Policy Detected is:

FeatureCount(%Feature Count%)
def FeatureCount(n):
if n >0:
    return False
else:
    return True

When there are features in the feature class (n>0), 'Policy detected' has a 'True' output and 'No Policy detected has a 'False' output - all correct!

When there are no features in the featureclass (n=0) then both 'Policy detected' and 'No Policy detected' show a False output instead of 'False and 'True'

Any ideas what's wrong?

.


texts: Locations of Geodetic Survey Benchmarks in Ontario?

texts: Does anyone know where I can find a list of geodetic benchmark locations for Ontario? (If a list even exists that is...)


texts: Merge to New Raster resulting a raster with cell value 0 to 255

texts:
I am using the mosaic to new raster tool to mosaic two rasters together, one 2-bit and one 4-bit, and the output raster has given all the previously nodata pixels the value 3. Some of the pixels in the output should be 3, their original and intended value, while the majority should be nodata or a much higher or lower value like 255 or -/+ 32,768 that appears when I use the mosaic tool on data with different extents.

Any idea what could be causing this? Should I be converting the 2-bit data to 4-bit to match the other raster?

I have used the tool on very similar data for the exact same process and not had an issue like this before.


texts: QGIS add new field for each value in exisiting field on a vector layer

texts:
I have a vector layer containing around 17000 polygons with a field named 'theme'. Each polygon of the layer has one of 30 unique values (e.g. buildings, roads etc.) assigned to this field. I now want to add a field to each polygon, the field should be named after the field 'theme' and contain the area of the polygon (e.g. a Polygon with the attribute 'buildings' in the field 'theme' gets a new field called 'buildings' containing it's area).
In Model Builder I can use the expression editor for the field name, my thinking was that I could get what I want with the expression: attribute(  $currentfeature ,'theme')
which does not add any new fields to the output at all (it does not generate any errors as well though).
Is this possible with Model Builder? Did I even use the correct expression or do I need to write some code to get this done?
The only way I found that works is creating 30 field calculators, one for each unique value and a lot of expression madness to only populate certain fields for certain values, since the Values might differ on a project basis this is really time consuming.


texts: Calculate the overlapping percentage two web feature service poygons

texts:
Is it possible to calculate the overlapping percentage of two polygons coming from web feature service?

I overlayed two WFS layers in my openlayers map.
The two layers have same data but collected at different timespan.So, even though they are same data, there is a change in the polygon area and shape.

I want to calculate the percentage of overlap of corresponding polygons in these two layers.

Is there any api to calculate for such cases?

I heard about Turf.js, intersect.js but couldn't get the right solution for my problem.


texts: Extract water bodies using Sentinel 2 - Google Earth Engine

texts:
I would like to extract water bodies (a lake to be specific) from Sentinel 2 images on Google Earth Engine. I know the Global Surface Water Explorer has pixel data on water-no water, but for my study area, the amount of usable images is very limited. Thus I wanted to branch out to different products.
To be specific I would like to create a time series with the extent of the lake. By calculating the area of the connected pixels for every month using the entire imagecollection. So far I know I have to use a function to go over every image, calculate the NDVI, mask all the pixels with a positive value and convert the raster to vector the calculate the area of the the neighbouring pixels with a negative value.
I'm having trouble putting all that into a working script, and would like to have some input on it.
// Create image collection of S-2 imagery
   var S2 = ee.ImageCollection('COPERNICUS/S2');
// Function to mask cloud from built-in quality band information on cloud
  var maskcloud1 = function(image) {
                   var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1));
  };
// Function to calculate and add an NDVI band
  var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B4']));
  };
// Function to mask cloud from built-in quality band information on cloud
  var maskcloud1 = function(image) {
                   var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1));
};
// Function to mask all non-water bodies
  var collection =  ee.ImageCollection.fromImages(
              addNDVI.map(function(image) {
  return image.updateMask(image.lt(0));
  }));
// Convert the zones of the thresholded to vectors.
  var vectors = collection.reduceToVectors({
  geometry: roi,
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'zone',
  reducer: ee.Reducer.mean()
  });

Anyone who has some insights?
Edit: clarification on the code. Right now I get an error when performing the mask function. With the code I tried to include Al the different aspects that I need. However, I'm not sure of their linked correctly.


texts: I would like to reverse Geocode from the infowindow

texts:
I am just learning the ropes with CartoDB and am curious to know if this is possible.

I have a CartoDB map with points and I am entering data into this map through the CartoDB editor.

I would like to be able to Google reverse geocode the street address of my point automatically.

Ideally I would create a link on the infowindow labeled 'Find address' this would then run a script which would find the address of that point and INSERT it into that points specified column on my table.

I would like this to run entirely from the CartoDB Editor as this is where I update my data from.

Does anyone know if this is possible and if so where is the best place to start?

If this won't work are there any other reverse Geocode solutions out there?


texts: From Kernel Density to a polygon that covers 80 % of points

texts: I have a set of occurrences as points in a district. I used Kernel Density tool, to show the areas of various densities. I also want to create one polygon that covers X % of points. I've tried the Locate Region tool but dont feel like this is quite the right method, as its not clear to me how it works and seems to be more for finding locations at set distances to other 'suitable' locations. I'm working in ArcGIS 10.1.5.


texts: Python - Projecting different shapefiles in the same folder

texts:
I have a few shapefiles in one folder with different coordinate systems but would like to use one of them as a template dataset. Is there a way to use python to filter out the shapefiles already in the required projection in the same folder?

I have been looking at the ListFeatureClasses and somehow using a wildcard, but not sure if that's the right way to go.


texts: How to find and replace all null values in multiple columns?

texts:
To find and replace all null values in a specific column I am using:

case when "myColumn" is NULL then '99999999' else "myColumn" end

But how can I use the QGIS calculator to find and replace all null values in the multiple columns?


texts: How to create a shapefile of circles from geo coordinates and radius in meters and project it correctly?

texts:
I have latitude and longitude coordinates for India from Google Maps in degrees and radius in meters as displayed below:

I want to create a shapefile of circles with given radius. I know about Analysis - Buffer tool, my problem is mostly about projection because the circles look very large, though the radius is half kilometer only. I don't know which projection to use to display the BUFF_DIST in meters and correctly display my data.
I tried Indian_1975_UTM_Zone_48N and the circle covers the whole India. I guess the problem is whether it assumes the radius in meters or km.
How can I project correctly Indian coordinates? I am using ArcMap 10.4.1


texts: Cannot plot OSM data retrieved from R package 'osmdata' with leaflet directly

texts:
I am using the R package osm_data to retrieve some polygons. If I want to plot them using leaflet, it does only work if I export the polygons to a shapefile and re-import it, but plotting it directly does not work.

Plotting using the provided plot variant does work, however.

I have played around with the format that osm_data returns, but to no success. Does anyone have an idea for the reason?

library(osmdata)
library(sf)
library(leaflet)
library(tidyverse)

result <- opq(bbox = 'Lisbon, Portugal') %>%
  add_osm_feature(
    key = 'amenity',
    value = 'university'
  ) %>%
  osmdata_sf()

universities <- result$osm_polygons %>%
  select(osm_id, name) %>%
  drop_na(name)

head(universities)

universities %>% select(name) %>% plot

leaflet() %>%
  addTiles %>%
  addPolygons(data = universities)

universities %>%
  write_sf('universities.shp', driver = 'ESRI Shapefile')

universities_from_shapefile <- read_sf('universities.shp')

head(universities_from_shapefile)

leaflet() %>%
  addTiles %>%
  addPolygons(data = universities_from_shapefile, label = universities_from_shapefile$name)

Output of head(universities)

> head(universities)
Simple feature collection with 6 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -9.209042 ymin: 38.65691 xmax: -9.136679 ymax: 38.74995
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
           osm_id                                                              name                       geometry
28681498 28681498                                 Faculdade de Medicina Veterinária POLYGON ((-9.195866 38.7140...
52430141 52430141 Faculdade de Ciencias e Tecnologia da Universidade Nova de Lisboa POLYGON ((-9.203055 38.6594...
71571983 71571983                           Universidade Católica Portuguesa Lisboa POLYGON ((-9.166331 38.7480...
93763775 93763775                                        Instituto Superior Técnico POLYGON ((-9.140742 38.7377...
95417328 95417328                   Instituto Superior de Ciências da Administração POLYGON ((-9.14661 38.70841...
97822065 97822065                  Universidade Nova, Faculdade de Ciências Médicas POLYGON ((-9.139408 38.7189...

output of head(universities_from_shapefile)

> head(universities_from_shapefile)
Simple feature collection with 6 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -9.209042 ymin: 38.65691 xmax: -9.136679 ymax: 38.74995
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

# A tibble: 6 x 3
  osm_id   name                                                                                geometry
  <chr>    <chr>                                                                                <POLYGON [°]>
1 28681498 Faculdade de Medicina Veterinária           ((-9.195866 38.71404, -9.195883 38.71407, -9.195884 38.71412, -9.1…
2 52430141 Faculdade de Ciencias e Tecnologia da Univ… ((-9.203055 38.65946, -9.203157 38.659, -9.203167 38.65862, -9.203…
3 71571983 Universidade Católica Portuguesa Lisboa     ((-9.166331 38.74801, -9.167548 38.74971, -9.167599 38.74991, -9.1…
4 93763775 Instituto Superior Técnico                  ((-9.140742 38.73776, -9.140703 38.73779, -9.140682 38.73782, -9.1…
5 95417328 Instituto Superior de Ciências da Administ… ((-9.14661 38.70841, -9.146233 38.7083, -9.146319 38.70812, -9.146…
6 97822065 Universidade Nova, Faculdade de Ciências M… ((-9.139408 38.71897, -9.139442 38.71903, -9.139504 38.71916, -9.1…


texts: ArcMap Merging multiple GeoTagged photo points

texts: I have a data set of Geo tagged photos that I imported into ArcMap using the GeoTagged Photos to points tool. Some of these photos were taken at the same location but because the Drone moved around a little bit with the wind each photo is represented with a new point that is only a few feet apart. How can I merge these photos so they are represented with a single point on the map?


texts: Can't get the correct geometry values after MySQL querying

texts:
I have large points data and also a small bounding boxes. I want all the points present within this bounding boxes.

With the extent of the bbox , I want to get all the points within this bbox.
With the help of the below code I was able to query and also write it to a .csv file.
The geometry column has non-geometric type variables when comparing the geometry column of the Queried_values and Original_values.

Required Output : get correct geometry values in the SQL queried data

This is my code below :

conn = MySQLdb.connect(host="", port=, user="", passwd="", db="db1", charset="")
cur = conn.cursor()
cur.execute("""SELECT * FROM db1.table1 where (latitude >=4.139916666666666 and latitude <=52.96272222222222) and (longitude >=5.000083333333333 and longitude <=51.999944444444445)""")
set1 = cur.fetchall()
with open("point.csv", 'w', newline='') as f_handle:
    writer = csv.writer(f_handle)
    header = ['latitude','longitude','geometry']
    writer.writerow(header)
    for row in set1:
        writer.writerow(row)

My Output:

Queried_values : This is the output I get after SQL querying.

Original_values : But I would like get the geometry as in this original file.


texts: Tracing or snapping to polygon border for measure tool

texts:
I am trying to do a simple measure of the length of seperate sides of a polygon depending on which other polygons it borders (in the image: the dark blue borders of the green polygon to the white area and the red border to the light blue area).

While I'd be keen to hear what tools would be useful for this exercise I have the problem that the polygons come from different sources and do not border exactly (gap) or overlap.
I am thus trying to measure the length using the 'Measure' tool in ArcGIS 10.5. This works fine for straight lines but the snapping of the measure tool to curved lines (holding down CTRL) is fairly imprecise. Is there any method to trace outlines similar to the 'Trace' tool or which other method could I use?


texts: How to connect a PyQGIS plugin with Postgres?

texts:
I have some questions for you:

I created a python-plugin in QGIS,

then I created a Qt gui fom my plugin

and at last I created a Postgres DB

now I need to make a connection between Qt-gui and Postgres-DB

I wrote this

 from PyQt4.QtSql import *
 db = QSqlDatabase.addDatabase("QPSQL");
 db.setHostName("localhost");
 db.setDatabaseName("siti");
 db.setUserName("postgres");
 db.setPassword("bina");
 bool ok = db.open();

but it doesn't works, I can't understand what I do wrong...

I wish someone could help me...

thank you for all!
lisa


texts: Raster Exporting or Clipping in ArcGIS Pro?

texts:
I don't do a lot of work with Rasters.

Basically I get a nice clean Raster when I use the "Topo To Raster" tool in 3D Analyst from a shapefile of contours.

This is a Rectangle (extent of my contours) so I need to clip it based on actual perimeter of my contours.  But no matter if I Clip using the smaller boundary (minimum boundary tool) or Export

using that new smaller boundary, I always get a lower quality Raster that has some pix-elation going on.

How do I create the same quality raster as the first, larger raster?


texts: Polygonize each object in layer individually (QGIS)

texts: Running the polygonize algorithm from the processing toolbox doesn't seem to work if you have a layer with several line objects that you want to polygonize individually. Think two "polygons-to-be" made of self-intersecting lines that also intersect each other. The algorithm takes the intersection of these two polygons as another polygon, creating three polygons in total although there are only two line objects in the layer. Is there a way to polygonize each object in a layer individually?


texts:
I tried to use the Raster (Extration) tool "clipper" and I always get this error log:
What is my problem?

Fehler bei der Ausführung von Python-Code:

Traceback (most recent call last):
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\GdalTools.py", line 342, in doClipper
    self.runToolDialog( d )
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\GdalTools.py", line 370, in runToolDialog
    dlg.show_()
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\tools\doClipper.py", line 47, in show_
    BasePluginWidget.show_(self)
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\tools\widgetPluginBase.py", line 44, in show_
    self.onLayersChanged()
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\tools\doClipper.py", line 71, in onLayersChanged
    self.inSelector.setLayers( Utils.LayerRegistry.instance().getRasterLayers() )
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\tools\GdalTools_utils.py", line 161, in getRasterLayers
    return filter( self.isRaster, LayerRegistry.layers )
  File "C:/Dokumente und Einstellungen/user/.qgis//python/plugins\GdalTools\tools\GdalTools_utils.py", line 156, in isRaster
    if layer.usesProvider() and layer.providerKey() != 'gdal':
AttributeError: 'QgsRasterLayer' object has no attribute 'usesProvider'

Python-Version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS-Version:
1.9.0-Master Master, 298f0c0


texts:
CartoDB helpfully shows a daily hit count for each map, but I can't find out how end of day is defined.  I seem to get lots of hits "overnight" in UTC, and I swear the hit count for two days ago just changed. What is the time zone?

Also, does just viewing a page with a map embedded count as a hit (eg by a robot), or does the visitor have to click into the map?


texts: QGIS: Error when using Undo/Redo on layers loaded from PostgreSQL database

texts:
I am working with layers loaded from PostgreSQL(11.7)/PostGIS database to QGIS 3.10.7. When I try to undo/redo changes made on the layer it fails, sometimes during the undo and basically every time when doing redo (I didn't notice any pattern in it). Apparently, there is some problem with the cursor in the PostgreSQL database, but I don't know why.
Fetching from cursor qgis_398 failed Database error: ERROR: cursor "qgis_398" does not exist

I thought that QGIS is able to handle in the background all the functionality for undo/redo on database layers (in newer versions at least).
UPDATE: It looks like when the layers are in editing mode, changes (creating new features at least) are happening on the QGIS side and are not propagated to the database. When I do "redo" QGIS is trying to obtain the data from database via cursor, but the cursor was not set beforehand. Probably because the changes were not pushed to database?
UPDATE 2: Over time I tried several versions of PostgreSQL (9, 11, 12) and QGIS (3.10, 3.14), the same problem persisted. After studying logs made by the database server I have found, that the QGIS is saving the state of the layers as a SAVEPOINT in database when I start editing. During the edits, QGIS is creating various cursors (with ST_AsBinary function). I am guessing that by doing this, the QGIS gets the actual data from the database. The problem occurs when the QGIS calls for the ROLLBACK (when I do undo in the QGIS). All the changes that was done after declaring the SAVEPOINT are lost, but QGIS then tries to FETCH the CURSOR that is lost thank to this. This causes the error and the transactio is aborted.
"DECLARE CURSOR",2020-11-24 15:51:20 CET,15/18,126847,LOG,00000,"duration: 0.473 ms  statement: DECLARE qgis_305 BINARY CURSOR WITH HOLD FOR SELECT st_asbinary(""geometry"",'NDR'),""id"",""fk_MK"",""popisky_z""::text,""popisky_y""::text,""poznamka""::text,""typ_useku""::text,""umisteni""::text,""provoz""::text,""popisky_x""::text,array_out(""par_cis_ve"")::text,array_out(""par_cis_hl"")::text,""KU""::text,""stav_MK""::text,""cislo_kom""::text,""cislo_us""::text,""povrch""::text,""plocha""::text,""delka""::text,""sirka""::text FROM ""MK"".""usek"" WHERE ""geometry"" && st_makeenvelope(-533970.59677669638767838,-1124501.15351219591684639,-532212.18808578583411872,-1123600.86123963561840355,5514)",,,,,,,,,"QGIS"
"ROLLBACK",2020-11-24 15:51:20 CET,15/18,126847,LOG,00000,"duration: 0.065 ms  statement: ROLLBACK TO SAVEPOINT ""qgis779d8f397d0d4932aa0a""",,,,,,,,,"QGIS"
"FETCH",2020-11-24 15:51:20 CET,15/18,126847,ERROR,34000,"cursor ""qgis_305"" does not exist",,,,,,"FETCH FORWARD 2000 FROM qgis_305",,,"QGIS"
"CLOSE CURSOR",2020-11-24 15:51:20 CET,15/18,126847,ERROR,25P02,"current transaction is aborted, commands ignored until end of transaction block",,,,,,"CLOSE qgis_305",,,"QGIS"


texts: Difference between GIS mapping and Photogrammetric mapping

texts: I can't seem to understand both type of mapping. I am planning to make a drone for mapping and I am confused about picking one. My purpose is to take photos of specific area and do 3D modeling of that area into the software.


texts: how do i assign values to shapes to be displayed in my map

texts:
I apologize for the simplistic nature of this question, I'm VERY new to QGIS. I'm trying to use QGIS to display political information. I have a map of my state, with each town as it's own shape. I have an ID for each town/shape. What I want to do is assign a value to each town/shape, and have a scaled display of these values by town/shape.

Can anyone help me with the process for assigning values to each town/shape?

Thanks!!!


texts: Merging rows based on duplicate field values using ArcGIS Desktop?

texts:
How do I merge rows with the same values in one of my fields in my attribute table.

I cannot do this manually due to the extent of the duplicates.


texts: How to activate and deactivate popup feature in OpenLayers2?

texts:
I am a beginner in OpenLayers and I want to activate and deactivate a popup feature, which has to show lat,long of the click, which is initiating the popup [no markers].

I am currently using OpenLayers.Control.WMSGetFeatureInfo and a variable to switch between getinfo and popup. But it has some issues, whenever I am clicking always it is making calls to WMS server, so I want separate controls, which I can enable and disable as I need...

Any ideas welcome...


texts: Using processing toolbox conditionally in QGIS?

texts:
I want to be able to run the 'split with lines' function, but only on polygons where the area > determined value.

Is it possible to do this simply, ie running a process from the processing toolbox only when a condition is met in QGIS?


texts:
I'm working in QGIS 2.18.3.

I have point layer which represents GPS points along a transect line. I would like to randomly select 5 points along this transect, i.e. from this layer.

Somehow I expected this to be straight forward, but I seem to miss something here.

I used Research Tools --> Random selection within subset, selected the transect as input layer, Latitude as ID Field (Latitude = Real), Number of selected features (for the count instead of percentage) and entered 5 as number of selected features.

After pressing "Run" the process runs without problems as the process bar proceeds to 100%, then closes and I receive no error message.

However, nothing has been selected, i.e. no 5 points are highlighted in the attribute table (header still shows selected=0).

Am I using the function in a wrong way or do I expect the wrong output?


texts: Getting ArcGIS to include leading zeros in text field?

texts:
I am entering characters in a text field that is essentially a 5 number ID. I need leading zeros to be included in the value, but every time I enter or calculate the field, the leading zeros are removed(even though the field properties say it is a 5 character text field).

How do I force the field to include leading zeros?


texts: Python bindings for PCL (Point Cloud Library)

texts:
I am curently searching for some Python 3.x bindings for PCL on Ubuntu (16.04).

The ones I found are:

pclpy: pip cannot find it.
python-pcl: I have tried to install it but I always end up on errors with missing files or unmet dependencies or version issues... What is strange is that it is the only one listed on the PCL main page: http://pointclouds.org/news/tags/python !
PyPCL: I can install it with pip but it's not documented at all.

Which one is the most used, the most efficient and up-to-date to process point clouds from Python?
Did I miss some other package(s)?

Edit:

to answer properly the 1st comment below:

I tried to install python-pcl:

1st error along with many many warnings about deprecated functions:

pcl/_pcl_180.cpp: In function ‘int __pyx_pf_3pcl_4_pcl_18ConditionalRemoval___cinit__(__pyx_obj_3pcl_4_pcl_ConditionalRemoval*, __pyx_obj_3pcl_4_pcl_ConditionAnd*)’:
pcl/_pcl_180.cpp:48508:153: error: no matching function for call to ‘pcl::ConditionalRemoval<pcl::PointXYZ>::ConditionalRemoval(__pyx_t_3pcl_15pcl_filters_180_ConditionBasePtr_t, NULL)’
__pyx_v_self->me = new __pyx_t_3pcl_15pcl_filters_180_ConditionalRemoval_t(((__pyx_t_3pcl_15pcl_filters_180_ConditionBasePtr_t)__pyx_v_cond->me), NULL);

My PCL version is 1.8.1.


texts: QGIS 3D Globe Plugin Elevation Map Projection

texts:
I'm trying to get a map setup using the QGIS 'Globe Plugin'. I've been trying to use the UK Environment Agency LIDAR data as an elevation map. After a bit of puzzling this out I've found that the Elevation Map needs to be in WGS 84 projection.

Is it just the case that you need to reproject the elevation map to be WGS 84 before using it, or is there a way to use other projections?

I hope that explanation isn't too confusing, but let me know if you need clarification on anything.


texts: Reprojection doesn't change the coordinates

texts:
Following the StackOverflow guidance, I have reprojected raster layers into WGS84/UTM 36N using the Raster-> Warp -> Reproject tool. According to the Raster Properties, all the rasters are now in the same CRS; however, the coordinate system is different and the rasters do not overlap (even though "On The Fly Projection" is on). Here are two screenshots of the problematic rasters with the different coordinates:

Based on comments from other threads, I understand that the units for the UTM projections is meters; however, the re-projection doesn't seem to properly changing the coordinate units. One remains in longitude/latitute while the other is in decimal degrees. Based on advice online, I have tried starting from scratch to ensure that the Map Project projection is correct, but nothing appears to work. Does anyone have any ideas? I am using QGIS 2.18.


texts:
Here is the plan. I need to delete every metro's least populated census tracts (in terms of population density). To what extent? As to preserve 98% of the metro's total population. I don't know how to program. Is there any existing tools I can use?

Thanks everybody!!

S.


texts:
I'm using OpenLayers 3 to interact with some maps. I first declare my map:

map = new ol.Map({
            target: 'map',
            layers: [
              new ol.layer.Tile({
                  source: new ol.source.OSM()
              })
                ],
            view: new ol.View({
                center: [0, 0],
                zoom: 2
            })
        });

I've an event that triggers an action, to change my map's view center. This way, (my coordinates are in EPSG:4326 format):

function CenterMap(lat, long) {
     console.log("Lat: " + lat + " Long: " + long);
     map.setView(new ol.View({
            center: ol.proj.transform([lat, long], 'EPSG:3857', 'EPSG:4326'),
            zoom: 5
     }));
}

When the function runs I get this at the explorer console:

Lat: 9.0412851667 Long: -79.5658145000

But the maps goes to [0,0], does anyone knows why this happen?


texts: What is it called a kind of Rubberband effect on screen when user moves the Feature participating in geometric network in ArcMap?

texts:
I am using Visual Studio 2012, C#, .Net Framework 3.5, ArcMap 10.2.1, ArcFM 10.2.1b.

I am working on the requirement as follows:
when user is in editing state in arcmap 'This Custom Tool' should be enabled. When user clicks the 'This Custom Tool' it will show a 'Custom Form' having Tab control with 2 tab pages. each tab pages will have Tree View and Previous|Next or Finish button to move UI forward/backward. The tab page1 is called Selected Features by default form loads with tab page 1 set to visible and focused. Then user will be able to draw a box on current focus map and retrieve all the point features which are from the editable feature layers and puts them on TreeView as Feature Layer->Feature Display Name (ArcFM ID8 List Item Display Name) on first tab page tree view.
On First tab page TreeView user will be able to check or uncheck the individual features or feature layers to move them to next tab page tree view this is just a filter process to avoid showing unwanted features in second tab page tree view.
When user clicks the Next button on tab page 1 it clones all the checked nodes and add it to second tab page tree view so that user can start editing the geometry for checked features.
Depending on current focus tab page (1 or 2) 'This Custom Tool' will either allow user to draw a box to select features from map or on click it will move the selected feature from the 2nd tabl page tree view.

I am able to set the new geometry and move the selected feature to new location without any issue but user wants to see the rubberband effect while moving the mouse on the map if feature is participating in geometric network, similar to what ArcMap edit tool does when moving the feature participating in geometric network it shows all the connected features movement on mouse move.

Can anyone point me out what is that rubberband effect called and how do i implement it in 'This custom tool'.

I did googled it with several different relevant phrases but nothing popped out which will help me look into the API to implement this effect in custom tool.


texts: Selecting by attributes on features sharing common value

texts:
I have a point vector with features sharing non-unique "ID"s. Each feature also has a "distance" field. I want to select among each set of features sharing the same "ID" only the feature with the minimum value.
Example:

row
ID
Distance

1
One
3.4

2
One
1.6

3
Two
5.4

4
Two
6.8

From the above, it should select rows 2 and 3.
I have seen similar questions that ask for the minimum value for each "ID",  but I want to select the whole feature (because I need the coordinates of the feature).
I also know that this problem is trivially handled with VLOOKUP in Excel, but I want to streamline my workflow and so I am looking for a solution within QGIS.


texts: adding points to a line shapefile at specified distance in python

texts:
I want to add create nodes\points at an interval on a line shapefile let say at every 100 meters and save those points as shapefile.

I came across a plugin for Qgis called Qchainage, this exactly performs the task i want to do, but I am not sure how to adapt this plugin into my python script.


texts: Is there any file that contains already made shield images that are frequently used and comes with the installation of TileMill, or is there any page I can download it from. Link would be appreiciated. Specifically, I need some already made images of shields for highway labels and house numbers.


texts: Preventing ArcGIS crashes when exporting Data Driven Pages multiple times?

texts:
I am using ArcGIS 10.3 for Desktop.

When I export a data driven MXD using a strip index (i.e. 3 pages into a single pdf) it works fine the first time I export it.

However, if I try and export the same MXD again it crashes and I have to close and re-open the MXD.

This only occurs when data driven pages are used.

Does anyone have a solution to this, or is this a bug within the software?


texts: Problem with font size depending on screen resolution - QT

texts:
I developed a QGIS plugin by designing the interface with QT Designer.

The inferface appears correctly on my computer (left image).

However, on some computers, I encounter size problems, especially font size (right image).

I noticed that the problem was encountered when in the windows settings "scaling and layout", the variable
"change the size of text, applications and other items" was set to 125% instead of 100% on my computer.

I have read a few things about the use of layout and the declaration of the environment variable QT_AUTO_SCREEN_SCALE_FACTOR =1.

Is there a trick to avoid this kind of problem?


texts:
I want to import shapefile into PostGIS. There is one database postgres in it but no table added yet.

I am using this command:

shp2pgsql -I -s 4326 districts.shp districts | psql -d postgres -U postgres

and it's giving me some error:

  Unable to convert data value to UTF-8 (iconv reports "Invalid or incomplete multibyte or wide character"). Current encoding is "UTF-8". Try "LATIN1" (Western European), or one of the values described at postgresql.org/docs/current/static/multibyte.html

Is my command wrong ??


texts:
In GeoServer, how can I enforce user to lock out after certain attempts?
I couldn't find way to enforce number of attempts when I create new user and set the password policy.


texts: How to loop through all columns in a table/layer getting all data for a row/feature in MapBasic?

texts:
I'm trying to extend this GNU-licensed plugin MI2SVG, which converts MapInfo maps to SVG, so that it writes associated table data into the SVG it creates, so that browser based JS web applications can read map data from MapInfo-created interactive maps.

It's going well except I can't find the MapBasic way to loop through all columns for one row in a layer, in a way that gives the index (column name or name) in each iteration alongside the data.

The only MapBasic loop examples I can find are for looping through table rows using the very specific Do While Not EOT(), which seems to be always used for rows and doesn't appear to give each iteration's index (or if it does, I can't see where).

So for example, here's code that loops through rows in a layer table i.e. features in a map and prints the contents of column 1 into a string:

sLayer = LayerInfo (nMap, nLayer, LAYER_INFO_NAME)

'Initialize alias variables
aObj = sLayer + ".obj"
aRow = sLayer + ".rowid"
aCol = sLayer + ".col1"

'Convert each feature in the layer to SVG
Fetch First From sLayer
Do While Not EOT(sLayer)
    objA = aObj
    nRow = aRow
    sCol = aCol
    'some code that prints SVG like "<path name="""&sCol&...etc

How would you loop through all the columns, getting index and values for each? My ideal output would be, instead of name="""&sCol&""" a string like this: data-col1="[col1data]" data-col2="[col2data]" data-col3="[col3data]" for tables with any number of columns, printing the index (name or number) and value concatenated with some other text for each column.

I'm guessing it would be something like this (but with ???s replaced with actual things I can't find in the MapBasic reference and User Guide):

Fetch First From sLayer
Do While Not EOT(sLayer)
    objRow = aObj
    sOutput = ""

    Fetch First From objRow
    Do While Not ???(objRow)
       sColumnName = ????Info ( someArgs?, ????_????_NAME )
       nColumnIndex = ?increment???

       'Not sure what type for this variable since data could be string or numeric?
       ?ColumnData =  objRow.col(nColumnIndex)

       outputString = outputString & " data-"&sColumnName&"="""&?ColumnData&"""

...or maybe instead of Do While Not ???(objRow) a for..next loop :

nColumnCount = objRow.???NumberOfColumns??? 'or' TableInfo( nLayer??, TAB_INFO_NCOLS )
for nCounter 1 To nColumnCount
   sColumnName = ????Info ( someArgs?, ????_????_NAME )
   ?ColumnData = objRow.col(nColumnIndex)
   'some code that prints SVG like "<path name="""&sCol&...etc
next

...but despite wading through the MapBasic reference, user guide and trying to fill in the gaps with http://testdrive.mapinfo.com/ I can't find what to use to fill those blanks.


texts: Faster st_is_within_distance for points to lines

texts:
I am looking for a fast way to run a st_is_within_distance (or equivalent) in R, to find all the lines (y) within a given distance of each point (x). (This would be an equivalent of ArcGIS's "neartable", for which we do not seem to have an easy equivalent in R's GIS packages, although there are packages for points-to-points, see e.g. Is there an R equivalent of NEAR in ESRI ArcGIS?) The R help for ?st_is_within_distance contains this note:

For most predicates, a spatial index is built on argument x; see http://r-spatial.org/r/2017/06/22/spatial-index.html. Specifically, st_intersects, st_disjoint, st_touches st_crosses, st_within, st_contains, st_contains_properly, st_overlaps, st_equals, st_covers and st_covered_by all build spatial indexes for more efficient geometry calculations. st_relate, st_equals_exact, and st_is_within_distance do not.

So my interpretation is that st_is_within_distance is going to be no faster than calculating all possible pairs of distances and taking the ones within the given cutoff. From my testing, I find this is true.
So I had thought, as a way to improve on the speed of st_is_within_distance, to use an approach based on a spatial buffer: buffer the lines out to the distance I want (15km in my case), then use st_within to see if the points fall into the buffers (and pick which buffers). Conceptually this is similar to what st_is_within_distance is doing, although this gives slightly different answers, probably because the buffer is only an approximation to a "true buffer". Worrying about this turns out to be not important though, because, much to my surprise -- and despite using a spatial index -- this approach is SLOWER than calculating all possible distances! (about 1.7 times slower in my example)
Questions:
(1) Is there a way (in theory and/or in R implementation) to use a spatial index to obtain a speed improvement for finding all lines within a given distance of a point? Or is calculating all possible distances the only sure way to accomplish this? (Note that in sf, st_nearest_feature uses a spatial index and is very fast, but that only obtains the closest line to the point, not all lines within a given radius of the point. This suggests to me that a spatial index could perhaps speed up the st_is_within_distance call.)
(2) Why, despite using a spatial index, is the st_within approach slower than calculating all possible distances? I am looking for any intuition to help understand why I was wrong.


texts: How can I crop an OpenStreetMap and export it using QGIS?

texts:
I want to crop an OpenStreetMap (a specific area, for example a town). I want to save and export this cropped area using QGIS (format could be .osm or .shp).

I have searched everywhere but I could not find a valid answer. Is this not possible? One can crop an OpenStreetMap directly but it is limited by the size of cropping area. I assume this is a fairly common practice. What are the alternatives available?

I am using QGIS 3.10 over Ubuntu 18.


texts: Setting linear units for Buffer distance from arcpy.GetParameterAsText()?

texts:
I am trying to run a buffer tool as part of script that I will both use as a standalone script and import into ArcToolbox. I'm just wondering, if I use GetParameterAsText for the buffer distance, how do I set the linear units? E.g. I want it preset so they only have to type 300 and the script takes that to be meters? Here is what I have so far:

# Buffer the Thames feature class to create the Flood Risk Area
in_features = "Thames"
out_feature_class = "Flood_R_A"
# Run the Buffer Tool
arcpy.Buffer_analysis (in_features, out_feature_class, arcpy.GetParameterAsText(5),    "FULL", "ROUND", "ALL", "")


texts: The number of residential addresses/houses per zip code for Toronto, Canada

texts: Would it be possible to obtain the number of residential address (or total number of houses) at the zip/postal code level (three digit level) for Toronto, Canada? Where might I find such information.


texts: Error in rasterToPolygon trying to get several coverage layers

texts:
I have a very strange problem:

I am trying to get a raster with coverage from 0 to 1 in each cell for types of vegetation, I started with a raster of vegetation type, but I needed a raster with 12 times less resolution, here's my code

library(rgdal)
library(raster)

#Load raster
birds<- readGDAL("C:/Users/corcoranbarriosd/Documents/new_bats/Rnew_bats/Birds/2007_Raster.img")
birds <- raster(birds)

#transform raster to a resolution 12 times lower
ra <- aggregate(birds, fact= 17, expand = FALSE)

#get a layer of coverage of agroforest, which comprises cells type 0, 11, 13 and 4
AgroFo <- rasterToPolygons(birds, fun=function(x){x == 0| x == 11 | x == 13 | x == 14})
AgroFo <-rasterize(AgroFo, ra, getCover = TRUE)

#tranform it to a scale of 0 to 1
AgroFo <- (AgroFo/100)

SecVec <- rasterToPolygons(birds, fun=function(x){x == 2 | x == 3| x == 8 })
SecVec <-rasterize(SecVec, ra, getCover = TRUE)
SecVec <- (SecVec/100)

#do it again for other types of veg
upland <- rasterToPolygons(birds, fun=function(x){x == 4 | x == 5 | x == 10 | x == 16})
upland <-rasterize(upland, ra, getCover = TRUE)
upland <- (upland/100)

wetland <- rasterToPolygons(birds, fun=function(x){x == 1 | x == 6 | x == 7 | x == 12})
wetland <-rasterize(wetland, ra, getCover = TRUE)
wetland <- (wetland/100)

It al goes well for AgroFo and for SecVec, but I get mistakes in Upland and getland in the step of rasterToPolygons.

This is the message i get:

Error in .Call(Polygons_c, srl, ID) :promise already under evaluation: recursive default argument reference or earlier problems?Error during wrapup: promise already under evaluation: recursive default argument reference or earlier problems?

I find it very weird, since the two first layers work fine, and the last two fail.

let me know of any clarification I can make to make my problem clearer


texts:
Just discovered QGIS and really impressed! I've not formal GIS training so rely on Google to figure stuff out.

I'm using the SHP files made available by the Canadian Government as part of their CANVEC+ set. The SHP files are made available by NTS map sheet.

The naming shows what of data is in the layer. I've a batch files that copies the layers I'm interested in out of a mapsheet dir.

I was trying to rename the Layer with a PyQGIS script. Examples on the internet got me to the point of being able to list the layers. If statements to match what I need. But have been unable to find if I can assign a new name value that way?

Is there a better way to do this, because I'm guessing this a pretty common problem?

Can I use the same script logic to assign a transparency (or maybe even a style)?


texts: How to import SVGs to QGIS without QGIS changing it?

texts:
I've been having problems with importing .svg-files to QGIS (LTR 3.10.5). This is how the image looks in Photoshop CC (then exporting as SVG).

And this is how the SVG translates into QGIS. Basically it loses the lines and the center circle gets black fill.

Here's the SVG code too.

<svg id="Ryhmä_2" data-name="Ryhmä 2" xmlns="http://www.w3.org/2000/svg" width="4000" height="2000" viewBox="0 0 4000 2000">
  <defs>
    <style>
      .cls-1, .cls-2 {
        fill: none;
        stroke: #000;
      }
      .cls-1 {
        stroke-width: 41.667px;
        fill-rule: evenodd;
      }
      .cls-2 {
        stroke-width: 20.833px;
      }
    </style>
  </defs>
  <path id="Muoto_1" data-name="Muoto 1" class="cls-1" d="M0,1000H1600"/>
  <path id="Muoto_2" data-name="Muoto 2" class="cls-1" d="M2400,1000H4000"/>
  <circle id="Soikio_1" data-name="Soikio 1" class="cls-2" cx="2000" cy="1000" r="400"/>
</svg>

How can I fix this problem? I've tried doing the SVG with Inkscape too, but the result in QGIS is the same.

Sorry for the finnish pictures, I hope you can get a hang of the problem!


texts: Merging Aqua + Terra MODIS image collections

texts:
I run into an error when merging 2 image collections.

Goal: I want to merge the (identical) 250m daily Terra and Aqua MODIS collections into one image collection to obtain one collection with 2 images per day:

// ------MODIS 250M DAILY ----------------------------------------------
var modis_250m_daily_terra =  modis_250m_daily_terra.filterBounds(aoi).filterDate(date_start,date_end);
var modis_250m_daily_aqua = modis_250m_daily_aqua.filterBounds(aoi).filterDate(date_start,date_end);`

Method to merge the collections:
>`// ------MERGE MODIS 250M DAILY ----------------------------------------------
var modis1 = modis_250m_daily_terra.merge(modis_250m_daily_aqua);
// sort by date
var modis1 = modis1.sort("system:time_start");
print('MODIS Surface Reflectance Daily 250m',modis1);
var size_modis1 = modis1.filterBounds(aoi).filterDate(date_start,date_end).toList(100000).length();
print('Total images in modis1 (250m daily) in current time frame: ', size_modis1);
print('Merged 250m MODIS',modis1);

This seems to work until in the following function, this code runs into an error "Invalid argument specified for ee.List(): sur_refl_b01", whereas this function used to work without the merge.

// ------CLOUDFUNCTION----------------------------------------------------
var cloudfunction = function(image){
  //filter the image on its high reflectance part
  var cloud = image.select('sur_refl_b01').rename('cloud');
  var cloud01 = cloud.gt(1000).lt(4000);
  image = image.updateMask(cloud01).addBands(cloud);
  var area = ee.Image.pixelArea();
  var cloudArea = cloud01.multiply(area).rename('cloudArea');
  image = image.addBands(cloudArea);
  var stats = cloudArea.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: aoi,
    scale: 250,
  });
  return image.set(stats);
};
var collection = modis1.map(cloudfunction);
print ('Clouds pixels', collection);

I tried a way to 'join' the collections instead, but this doesn´t give the semidiurnal image collection. Any ideas on how to either avoid the error, or to improve the merge that I applied?


texts: Mean Distance Between features on Two Feature Class

texts:
I am new to arcmap and need your help to resolve the below issue.
I have two feature layers about a town which was moved from one location to another. Now i need to calculate the mean distance between the new and old town.
I tried using Euclidean Distance but i understand the value returned are not correct. I also have a digital elevation model of the area.

Please refer the attached screenshot


texts: How to extract Metadata from File Geodatabase using Ogr2ogr?

texts:
The Ogr2Ogr doc says:

GetLayerDefinition a_layer_name” and “GetLayerMetadata a_layer_name” can be used as special SQL requests to get respectively the definition and metadata of a FileGDB table as XML content (only available in Geodatabases created with ArcGIS 10 or above)

I'm able to connect, since this works:
ogrinfo -so -q -nogeomtype /path/to/myfile.gdb

However when I try this:
ogr2ogr -dialect OpenFileGDB -sql "SELECT GetLayerMetadata(myfeatureclass)" /path/to/myfile.gdb

It says: FAILURE: no source datasource provided
With ArcGIS Pro, I can see the metadata for the featureclass.
Are there any examples for this somewhere?


texts: Are OSM ids unique over all Object types?

texts:
Are OSM IDs unique over all Object types?

Meaning, can a way-ID occur as a relation-ID or a node-ID as a way-ID?


texts: What is the difference $GPRMC and $GPRMCA?

texts:
What is the difference between $GPRMC and $GPRMCA strings?

I can find online documentation about $GPRMC strings, however I have a device that sends $GPRMCA strings.

$GPRMC,225446,A,4916.45,N,12311.12,W,000.5,054.7,191194,020.3,E*68
$GPRMCA,101959,A,0,N,0,E,0,0,20160503,,*1A

The most obvious difference is that the day field has a different formatter (YYYYMMDD in $GPRMCA instead of the DDMMYY FROM $GPRMC)


texts: Automating Raster Clipper for single raster but multiple masks using PyQGIS?

texts:
I have a large raster that I need to clip into many smaller tiles, however as there are nearly 200 tiles it will take too long to clip each tile individually and I need to automate the process.

The following python script will clip multiple rasters with a single extent box however, I need to clip a single raster with multiple masks.

Is it possible to simply adjust the parameters of this script to accomplish what I need?

# input variables

# where to save clipped rasters
outputDir = "/home/user/Desktop/clip"

# upper left point
UX, UY =  xxxxxx , xxxxxx

# lower right point
LX,  LY = xxxxxx , xxxxxx

# actual script
import os
lddLrs = qgis.utils.iface.legendInterface().layers()
for lyr in lddLrs:
 if (lyr.type()==QgsMapLayer.RasterLayer):
  os.system("gdal_translate -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr.source()) + " " + outputDir + "/" + str(lyr.name()) + ".tif")
  qgis.utils.iface.addRasterLayer(str(outputDir + "/" + lyr.name() + ".tif"), str(lyr.name() + "_clip.tif"))


texts: Georeferencing with shapefile instead of raster image in QGIS

texts:
I have a geographic file (be it shapefile or otherwise) which is stretched and dislocated into the wrong projection space.
So I'm wondering how to approximate its correct location by selecting points on this file which correspond to real-life locations in the given basemap, and use this knowledge to correctly geo-reference it in the way we would with most raster images?


texts: Subquery a list in select layer by attribute management using python

texts:
I am new to python and having trouble with my script.

I want only the records from the COMTRS field (of the pur_dbf_view) that match values from the ComtrsList to be selected; however, I am getting an "invalid expression" error message.  I think to solve my problem, I need to write the python list (ComtrsList) to a .dbf or a file format that can be later referenced in the Select Layer by Attribute management at the end of the script.  Is this right?  Is my syntax for the subquery in the Select layer by Attribute statement correct?

import arcpy, sys, os, traceback

# Set workspace
arcpy.env.workspace = "H:\\GIS_AIR\\Python\\PUR_tool"

# Set variables
plss = "Plss.shp"
comm = "Communities_sub.shp" # Input Community feature class
c_layer = "Communities_layer_sub" # Name of the Make Feature Layer output for    communities
p_layer = "Plss_layer" # Name of the Make Feature Layer output
plss_out = "Output_plss.shp" # Name of PLSS (selected) output feature.  Make sure this    has the .shp extension

search_distance = '0 MILES'

pur_dbf = "Glenn_pur.dbf"
pur_dbf_view = "Glenn_pur_view.dbf"

pur_dbf_sum = "Glenn_pur_sum.dbf"

comtrs_list_dbf = "Comtrs_list.dbf"
template = "TemplateA.dbf"

# Make feature layer for PLSS feature
arcpy.MakeFeatureLayer_management(plss, p_layer)

print "Completed make feature layer management PLSS"

# Make feature layer for Communities feature
arcpy.MakeFeatureLayer_management(comm, c_layer)
print "Completed make feature layer management Communities"

# Select features that are within 0 miles of comm fc
arcpy.SelectLayerByLocation_management(p_layer, "WITHIN_A_DISTANCE", c_layer,     search_distance, "NEW_SELECTION")

arcpy.CopyFeatures_management(p_layer, plss_out)

# Create new table with template
arcpy.CreateTable_management("H:\\GIS_AIR\\Python\\PUR_tool", comtrs_list_dbf, template)

ComtrsList = []
ComtrsCursor = arcpy.SearchCursor(plss_out)
for row in ComtrsCursor:
    comtrsValue = row.getValue("CO_MTRS")
    ComtrsList.append(comtrsValue)
print ComtrsList

###Write ComtrsList to comtrs_list_dbf (.dbf)###
arcpy.MakeTableView_management(pur_dbf, pur_dbf_view)
print "Created table view"

# Select layer by attribute on the pur_dbf_view
arcpy.SelectLayerByAttribute_management(pur_dbf_view, "NEW_SELECTION", "\"COMTRS\" IN (SELECT \"COMTRS\" FROM comtrs_list_dbf)")

arcpy.Statistics_analysis(pur_dbf_view, pur_dbf_sum, [["LBS_AI", "SUM"]], "CODE")


texts: Building a Python cron job/web API for QGIS

texts:
I have some tasks that I'd like QGIS to do every day, or when triggered by an API call.

I'm trying to understand the architecture.

I can run QGIS on a VM, then have a Python script accessing it and an external database to process any imagery that has come in since the last time the script ran, right?

Alternately, is there a straightforward way for me to set up a Python web server and then take data from POST calls and run it through QGIS?


texts: Dividing polygon into specific sizes using ArcGIS Desktop?

texts:
I have several thousand irregularly shaped polygons in a shapefile. I want to be able to split each polygon into three areas, and to specify what the size of those areas (they sum to the previous total area). It does not matter what the sub-polygon shape is, as this is for visualization purposes.

How can I go about doing this? Is there a standard algorithm that exists that I can use?

One approach I considered was to get all the points that make up the polygon. Then, I would join two randomly together using a straight line, split the polygon, then check if the area was within a satisfactory tolerance. If it was too small, I would change the point in one direction; if it was too large, I would change to a point in the opposite direction.


texts:
Given a line on the earth's surface, how do I plot a line perpendicular to it?

Apologies if this is a very simple question. I thought this would be a straightforward task, but it is proving counter-intuitive.

I start with the blue line in the figure below (see link -- I'm unable to upload figures yet). I found a perpendicular line by calculating the gradient of the blue line (m), then plotting another line (green) with gradient -1/m. When I plot the lines in Matlab (using 'plot' and 'axis equal'), they look perpendicular, as expected.

http://imgur.com/7qMkx

However, when I export these lines to Google Earth (using the KML Toolbox), they no longer look perpendicular (see link below; the shorter line is the blue line from the previous figure).

http://imgur.com/ncJQ7

I understand that strange things happen on curved surfaces, but I thought the lines should at least look perpendicular locally. I suspect this has something to do with the projection in Google Earth -- in particular, the fact that the grid cells appear to have roughly similar side lengths, yet the longitudinal edge has length = 1 degree, whereas the latitudinal edge has length = 0.5 degree.

So, in summary:

is my method of finding a perpendicular line valid on a curved surface? (i.e., plotting a line with gradient -1/m)
in the Google Earth image, do the perpendicular lines look as expected, or is something strange going on?

UPDATE:

To provide more context: I am looking at radar data taken from a plane. The multi-colored area is the 'swath', where observations have been recorded. The blue line I started with in the explanation above is parallel to the swath: this is the aircraft's flight line (the plane was moving in the roughly south-west direction). The radar looks in the direction orthogonal to the flight line, on the left. I am trying to draw a line perpendicular to the flight line; this should be the direction the radar is looking, and should cut the swath neatly. As you can see, this is not the case.


texts: GeoServer Basic Auth Login Failures (but everything works)

texts:
I've set up GeoServer to use basic auth (oauth/oidc is a can of worms I dont like opening). I am testing it using QGIS. Everything seems to work fine in QGIS however I am getting a phenomenal spam of login failures in the GeoServer console.

26 Oct 18:01:10 WARN [geoserver.security] - Failed login, user MyUser from LANIP, forwarded for REMOTEIP, stopped 30concurrent logins during authentication delay

I am using GeoServer behind a simple nginx reverse proxy hence the LAN/Remote IP.
While it does work,I wanted to ask anyway if anyone had any idea.


texts:
Using GeoServer 2.8.4, I have been successful at creating JPEG2000 image mosaics with OpenJDK 1.7.0. I use the gdal plugin and jp2k plugin and JPEG2000 images are loaded with Kakadu. I explicitly set it.geosolutions.imageio.plugins.jp2k.JP2KKakaduImageReaderSpi when creating the Image Mosaic.

As far as I know, the JPEG2000 coverage store is made available by the GDAL plugin which uses other plugins, in this case the JP2K plugin.

But now I need to switch from OpenJDK 1.7.0 to Java HotSpot JRE 1.8.0 and I get the following message when creating the layer:

HTTP/1.1 500 Internal Server Error
Illegal grid range [0 .. -33,793] for dimension 0. : Illegal grid range [0 .. -33,793] for dimension 0.
---- Debugging information ----
message             : Illegal grid range [0 .. -33,793] for dimension 0.
cause-exception     : java.lang.IllegalArgumentException
cause-message       : Illegal grid range [0 .. -33,793] for dimension 0.
class               : org.geotools.coverage.grid.GridGeometry2D
required-type       : org.geotools.coverage.grid.GridGeometry2D
converter-type      : org.geoserver.config.util.XStreamPersister$GridGeometry2DConverter
line number         : 48
class[1]            : org.geoserver.catalog.impl.CoverageInfoImpl
converter-type[1]   : org.geoserver.config.util.XStreamPersister$CoverageInfoConverter
version             : 2.8.4
-------------------------------

Logs don’t show any error prior to this one. And logs don’t show any relevant difference between OpenJDK and HotSpot JRE.

I should mention that GeoServer runs under WildFly 8.2.0 and, besides this issue, everything works like a charm with HotSpot JRE 1.8.0.


texts: Running QgsZonalStatistics in PyQGIS script?

texts:
I have experienced a very bizarre issue with my PyQgis code. I simply want to load a raster and a vector layer to do zonal statistics, the first part works perfectly but for the zonal statistic I have got an error with this command:

raster file = qgis.utils.iface.mapCanvas().layer(1).source()

(AttributeError: 'NoneType' object has no attribute 'source')

The strange thing is that if I rerun only the second part of the code (ie. my layers are already loaded) the extension .source() does not report any errors. This is still annoying since I am not able to write a loop with the full code and then use it for different layers.

Did anyone already experience this kind of issue?

Here is the code,

os.chdir("C:\Projet")
###          Part I
##           load vector layer
vbaselayer = QgsVectorLayer("county_save_as_lamhert_onlyUS.shp", "basemap", "ogr")
##           load and project the raster layer

os.system('gdalwarp -overwrite -t_srs EPSG:102004 -of GTiff NETCDF:"rhum.2m.2000.nc":rhum rhum_projection.tif')
humidity = "rhum_projection.tif"
fileInfo = QFileInfo(humidity)
rlayer = QgsRasterLayer(humidity,"hum2000")
rlayer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(rlayer)
QgsVectorFileWriter.writeAsVectorFormat(vbaselayer,"2000.shp","utf-8",None,"ESRI Shapefile")
year = 2000
vlayer = QgsVectorLayer("{}.shp".format(year), "county", "ogr")
vlayer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(vlayer)
###          Part II
##           zonal statistics
id = 1
while id < 365:
    vectorlayer=qgis.utils.iface.mapCanvas().layer(0)
    rasterfile = qgis.utils.iface.mapCanvas().layer(1).source()
    zonalstats = qgis.analysis.QgsZonalStatistics(vectorlayer,rasterfile,str(id),id)
    zonalstats.calculateStatistics(None)
    id += 1


texts: How to get Network Analyst data from script file back into Model Builder?

texts:
I am running a script which takes some input parameters and runs some Network Analyst processes. The script is imported into a toolbox and is setup in Model Builder with the necessary input and output parameters.

My problem is that I am struggling to understand how to get the resulting Solve operation back out of the script and into Model Builder for subsequent work.

Here is my script:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("Network")

# Script arguments
Network_Dataset = arcpy.GetParameterAsText(0) # Network dataset layer

Input_Locations = arcpy.GetParameterAsText(1) # Layer file

Output_Path = arcpy.GetParameterAsText(2) # Output path

# Local variables:
#Solve_Succeeded = Network_Analyst_Layer__3_
#Output_derived_feature_class = Network_Dataset

arcpy.AddMessage("Creating the Origin Destination layer...")

# Process: Make OD Cost Matrix Layer
arcpy.MakeODCostMatrixLayer_na(Network_Dataset, "odMatrix", "Length", "100", "", "", "ALLOW_UTURNS", "", "", "", "NO_LINES", "")

arcpy.AddMessage("Calculating Locations. This may take a while...")

# Process: Calculate Locations
arcpy.CalculateLocations_na(Input_Locations, Network_Dataset, "500 Meters", [["Streets","SHAPE"],["Streets_ND_Junctions","NONE"]] , "MATCH_TO_CLOSEST", "", "", "", "", "", "", "", "", "", "INCLUDE", "")

arcpy.AddMessage("Adding locations to Origin Destination layer...")

# Process: Add Locations to Origins Sub-Layer
arcpy.AddLocations_na("odMatrix", "Origins", Input_Locations, "", "500 Meters", "", "", "MATCH_TO_CLOSEST", "APPEND", "", "", "INCLUDE", "")

# Process: Add Locations to Destinations Sub-Layer
arcpy.AddLocations_na("odMatrix", "Destinations", Input_Locations, "", "500 Meters", "", "", "MATCH_TO_CLOSEST", "APPEND", "", "", "INCLUDE", "")

arcpy.AddMessage("Starting Network Analysis solve...")

# Process: Solve
arcpy.Solve_na("odMatrix", "SKIP", "TERMINATE", "")

arcpy.AddMessage("Network analysis completed.")

# How to get this out of script???

# Return as model output parameter
arcpy.SetParameterAsText(3, outputTable)


texts:
I am using some shapefiles in an agent-based model I am developing using MASON (a Java-based modeling framework). All the shapefiles I am using have a EPSG: 27700 OSGB 1936/British National Grid CRS.

I have just started to test my model using a different study area and the model is visualising the shapefiles with a different CRS to what appears in QGIS/ArcGIS.

The shapefiles are clipped from the same larger shapefiles that I have used for the other study areas and those display as expected using EPSG: 27700.

I've been told that:

  GeoMason is "blind" to the CRS. It simply imports the geometric
  character of each feature, disregarding projections and whatnot. There
  is no problem with this approach, you just have to make sure every
  data set you feed to a particular simulation refers to a single CRS.

  GeoMason exports rasters and vectors without setting CRS information,
  although it is tacitly the same CRS as the input data. A pasteriori
  you might hardcode the CRS yourself using a tool like GDAL/OGR (note
  that this is not possible for certain formats like ArcGRID).

I have no idea how to "hardcode" the CRS myself and have never used GDAL/OGR and only have limited Python experience.


texts:
I have a river network (lines) and a .shp with states (polygon). I would like to split the lines at the state boundaries i.e. split lines by polygon boundaries.
Is there a geoprocessing tool to do this in QGIS? I can't merge the lines by the way, since they are meaningfully split into tributaries/main river.


texts: Display frequently updating data in Cesium directly, or in a Geoserver layer?

texts:
I need to show a map layer with polygons showing the geographic extent "coverage" of an area our sensors' data points are coming in from. I need to update that once a minute or sometimes more frequently as the sensors move, so we know the current geographic areas the sensors are collecting from. We are serving basemap and static map layers from Geoserver and displaying them on a Cesium map. We need to update the view in Cesium regularly so we can monitor where the sensors are.

Should we attempt to use Geoserver to host the dynamically updating "coverage" layer and retrieve it periodically for display on a Cesium map, or push the data to the web UI and draw a layer directly in Cesium?


texts: Seamlessly Repeating a VectorLayer in OpenLayers

texts:
I've created a map in VectorLayers and done a number of things like allow selection, removal of features, adjusting the view, etc. however, one thing that has me stumped is that unlike all the examples I've seen, my map doesn't repeat seamlessly along the horizontal axis.
I have a version running in code sandbox here:
https://codesandbox.io/s/seamlessly-repeating-vectorlayer-xw9fh
Note that if you click and drag the map horizontally it doesn't repeat.
Is there a property in the geojson file I should change or is it something in my setup?
I've also included my main js below:
import "ol/ol.css";
import Map from "ol/Map";
import VectorLayer from "ol/layer/Vector";
import VectorSource from "ol/source/Vector";
import View from "ol/View";
import GeoJSON from "ol/format/GeoJSON";

import * as worldMap from "./world.json";

var geoJSON = new GeoJSON();
var vectorFeatures = geoJSON.readFeatures(worldMap.default);
var vectorSource = new VectorSource({
  features: vectorFeatures,
  format: geoJSON
});
var vectorLayer = new VectorLayer({
  source: vectorSource
});

var map = new Map({
  target: "map",
  view: new View({
    center: [0, 0],
    zoom: 18
  })
});

// dynamically added in so it can change during runtime
map.addLayer(vectorLayer);


texts:
I am working to download the MODIS data for every 5 days for one specific area.
Is there any ways to download it by once but not manually select every 5 days?
I have tried download from EO browser and getting bands 2,4,5 and to calculate from ArcGIS to get the EVI for my ROI.
Is there any possible way to download them in once?


texts: Automatically download Modis image using python script

texts:
I need download Modis imagery automatically from the website: https://earthdata.nasa.gov/.

I found this link: http://gis-techniques.blogspot.com/2010/06/getting-modis-image-automatically-from.html

but I can't run it. In that code there's the line:

logFileread=open(r"D:\MODIS\MOD09\lpdaac.txt",'r')

I don't know the structure of this file.I don't know how to create it.
If anyone knows the answer, please tell me.


texts: Zonal Statistics As A Table missing Features in Results?

texts:
I have a point shapefile that includes parcel centroid features across a county in Pennsylvania.

I am trying to get a sense of the percentage of forest land within half-mile and mile buffers around each of the centroids. Using NLCD land cover data, I've created a raster that indicates the presence of forest (== 1).

I've then applied Focal Statistics to that raster to create a raster that indicates the % of area within a 1 mile buffer of each cell that is in forest.

I've then tried to apply these raster values to the parcel centroid shapefile using Zonal Statistics by table.

Unfortunately, it hasn't worked. Out of 28000 centroids, it's only working for 25000 of the points.

Why might this be happening?


texts:
I uploaded a custom shapefile (converted to geoJSON) using the mapbox python SDK and visualize it in a Jupyter Notebook using the MapboxGL python module.

After I successfully upload the GeoJSON file as a tileset the minimum zoom level is set to 5 which makes visualizing this data at a global scale impossible. Mapbox point to this doc and suggests going the tippecanoe route. Is there another option that doesn't require to compile another package? Ideally I would like to overwrite these zoomlevel restrictions if necessary.

Upload script