texts: Import .gpx data to PostGIS using ogr2ogr

texts:
I want to import a few .gpx files to postgresql using ogr2ogr.
The data in each .gpx looks like this:
<trkpt lat="41.842607000" lon="12.490628000"><ele>40.5</ele><time>2010-05-31T14:12:11Z</time> </trkpt>

Firstly, I create a new table in postgresql
CREATE TABLE test (
           fid serial NOT NULL, --
           geom geometry(Point,4326),
           ele double precision,
           "time" timestamp with time zone,
           CONSTRAINT activities_pk PRIMARY KEY (fid)
);

Then when I run the ogr2ogr script:
ogr2ogr -update -append -f "PostgreSQL"  PG:"host="xx" user="xx" dbname="xx" password="xx"" 1st.gpx  -nln test

However, the script failed due to geometry type mismatch.

Warning 1: Geometry to be inserted is of type Multi Line String,
whereas the layer geometry type is Point. Insertion is likely to fail
ERROR 1: ERROR:  Geometry type (MultiLineString) does not match column
type (Point)
ERROR 1: INSERT command for new feature failed. ERROR:  Geometry type
(MultiLineString) does not match column type (Point)

Can anyone help with this?


texts: Creating all possible polyline routes connecting many point locations around polygon obstruction using ArcGIS Desktop?

texts:
Lets say I have 30 ports on the west of an island, and 30 on the east, and need to calculate all possible shipping distances between each and every port, with the routes having to steer around the geography of the island, which I have as a polygon.

I could digitise by hand but this would take an immense amount of time, so ideally I would like to get ArcGIS to calculate these routes for me by connecting all the points with polylines, but with the polygon as a barrier the polylines have to avoid rather than simply connecting in a straight line over the land (polygon).

As far as I can see in Network Analyst, if I had already had a network of potential routes, it would give me the shortest route option and take into account the polygon barrier. I could create a basic 'rough' network based on a vector grid, and get it to give me routes from this, but is there a more elegant solution? It doesn't seem like it should be an unusual task, but I can't find an existing solution?

Any ideas?


texts: Updating field aliases for ArcGIS map document layer using ArcPy?

texts:
I am on ArcGIS 10.4.1.

I have an .mxd file with a bunch of layers. The underlying database table has got new aliases for certain fields and I would like to get this reflected in the fields aliases for the corresponding layers.

If I remove the layer and add the feature class, I have to re-set the symbology + all other layer properties such as scale ranges, labels etc.
If I save the .lyr file first and re-apply it using arcpy.mapping.UpdateLayer, I will get back all the symbology, but also the field aliases which I don't want.
If I set the symbology_only=True, I will need to re-set all other properties manually.

I have seen the ArcObjects code and am able to update it using comtypes and C# code, however I wonder whether there is a way to update the layer's field aliases using arcpy.

I have seen that ArcMap X-Ray add-in has functionality of updating field aliases for the map document layers, great one. However, this should be done for every map document manually, thus it's impossible to script.


texts: Should I exclude coastal aerosol band of Landsat 8-OLI while computing PCA?

texts:
I am using Landsat-8-OLI Multispectral images to perform Principal Component Analysis. The aim is to detect noise and reduce dimensionality contained by the spectral bands. The resultant Principal components having new data dimensions with minimum noise will be used in biomass study which involves Multi Regression Analysis. The aim is to ascertain whether there is any potential relationship between biomass samples and PCA components.

The OLI sensor consist of 7 Multi spectral bands including ultra blue (B1) , wavelength spanned between 0.43 - 0.45 micrometer.

https://landsat.usgs.gov/what-are-band-designations-landsat-satellites

Should I exclude Band 1 - Ultra Blue (coastal/aerosol) and run PCA analysis on remaining Six bands (B, G, R, NIR, SWIR 1 & SWIR 2)?


texts: What is the proper syntax and usage for ArcGIS "in_memory" workspace?

texts:
I am unclear how best to use the in_memory workspace in ArcGIS.  The best resource I have found on the subject is from ESRI on Using in-memory workspace.  Unfortunately, the help section is not very informative.  Also, there is a post on this site that briefly touches on the subject here and here.

My specific questions:

What is the proper usage and syntax for using in_memory workspace in
ArcGIS/arcpy based scripts?
Is in_memory workspace the same as, for example, creating a layer
using arcpy.MakeFeatureLayer_management()?
Are there any standards such as deleting in_memory workspace at the end of
the script?


texts: Using large-scale maps/plans in Esri Story Maps

texts:
If there´s a possibility to adjust the scale of the basemap to plan-like maps.
Thing is, that i´ve got a map/plan of a greenhouse and i want to make a "guide" through this greenhouse to show some rare species of plants. But there is a restriction in Story Maps builder and i can´t zoom closer to the greenhouse, so it´s so small.

You can take a look at FLORA | Subtropický skleník which is the story map.


texts: Importing theodolite recordings into QGIS?

texts: I recently mapped an archaeological site using a theodolite (or dumpy level). Therefore, each point recorded has a distance measurement from the theodolite and a degrees off north. I'd like to import this information into QGIS. I know how to import data from Excel but how do I tell  QGIS that the one recording is distance and the other direction?


texts: Calculate area inside polygons based on slope

texts: I have some DEMs and have computed the slope in degrees for each one. I have also classified the slope raster into the following 5 classes: 0-3, 3-15, 15-35, 35-50, >50.  What I want to do is to calculate the area that has, for example, slope inside the 0-3 degrees class.


texts: Leaflet map altering page scroll behavior

texts:
When scrolling down a page and past a leaflet map, if the user's mouse cursor is inside of the map, then the scroll bar stops and the browser window becomes stuck on the map.

Is there a way to completely disable the scroll behavior on the map so that scrolling operates as normal?

I'm already using map.dragging.disable();


texts: dijkstra_sp runs without error in postgres query terminal, but triggers an error in Java

texts:
I am running this simple query

select id,gid from dijkstra_sp('roads',22,33)

This works perfectly fine and returns results in postgres query terminal. However running exactly the same query in java triggers the error:

  org.postgresql.util.PSQLException:
  ERROR: type "geoms" does not exist
  Where: compilation of PL/pgSQL function "dijkstra_sp" near line 6

Could anyone tell me how to solve this issue?

yours sincerely.


texts: How can I change the projection used by the QGIS Browser preview?

texts: QGIS Browser has no settings in the UI. Is there any file or so where I can define this?


texts: Creating points based on coordinates in MapBasic

texts: I am currently working on a MapBasic tool that generates start and end coordinates for a polyline layer. So far I have successfully populated 4 columns in a table (SX, EX, SY, EX) with coordinates of start and end nodes. The tool then exports the table as a csv. Now I want to open the exported table (I can do this) then create points based on these 4 columns in mapbasic. How would I go about doing this?


texts: Filling spaces between polygons using ArcGIS Pro?

texts:
I have a polygon with 4 attributes (1,2,3 or 4).  In the process of creating this polygon (merge/erase/union), there were gaps left between 1 and 2,  with some of the gaps having slivers of 3 or 4.

How do I snap polygon 1 to polygon 2,  without moving polygon 2, and simultaneously remove any slivers 3 and 4, using ArcGIS Pro?


texts:
I'm attempting to import the 45GB (compressed) planet-latest.osm.bz2 file downloaded from ftp://ftp.spline.de/pub/openstreetmap/planet/planet-latest.osm.bz2 but I receive parsing errors. As of now I have had two separate files fail in the same way - downloaded approximately 1 month apart.

The import command is as follows:

nohup osm2pgsql -l --bbox -180,-59,-26,85 -U gis -c -I -v --unlogged --drop -s data/planet-latest.osm.bz2 &

and the errors I receive are:

...
Entity: line -1123806515: parser error : Attribute user redefined
er="Schandlers" uid="51:48Z" version="1" changeset="29" user="OrcaDan" uid="411"
                                                                               ^
Entity: line -1123806515: parser error : Attribute uid redefined
er="Schandlers" uid="51:48Z" version="1" changeset="29" user="OrcaDan" uid="411"
                                                                               ^
data/planet-latest.osm.bz2 : failed to parse
Error occurred, cleaning up

This is a real pain as the file takes several days to download and the import can be running for several more days before the error is thrown.

I'm trying to understand what is causing the failure and what I can do to avoid it. My only theory is that the file is getting updated on the server as I am downloading it and I somehow get an invalid combination of two files. However I'm not really sure this is possible...

My question has three parts:

what is causing the error I'm seeing?
is there a faster way to validate a planet file than waiting for the import to fail?
can anyone point me to a recent and known-good file, if the file itself is the problem?

osm2pgsql version: osm2pgsql SVN version 0.82.0 (64bit id space), installed from Ubuntu 14.04 repos.


texts: What is the recommended format to import tables from Excel to arc map with more than 65 536 rows? CSV?

texts: What is the recommended format to import tables from Excel to Arc map with more than 65 536 rows? CSV?Txt?


texts: Interpolating lots of points with seemingly good spatial density

texts:
I'm using Arcmap 10 and I have scattered points (denser in an area). Each point represents a thickness value of a simulated avalanche deposit. I want to get a raster from this points which I will transform in polygons (raster-to-polygons).
Which is the best interpolation tool for this operation?

I tried them and I saw the "kriging" tool is the one which gives me the best "overlapping"* (i.e., the values of the points and the resulting raster coincide). Is this an appropriate way of establishing which tool is the one that suits better to my case?

Below, there is the image of the points superimposed on the raster obtained from them with the kriging interpolation.


texts: Can I avoid the offline tile cache limit (i.e., 6k) by self-hosting a Mapbox vector tiles?

texts:
BACKGROUND
I'm developing an off-line native Android mapping app. All map tiles (prefer vector, raster as fall back) must be cached locally on the device. Then the map will be overlaid with points of interest (POI), lines/shapes. The app will allow users to add their own POIs and draw lines/shapes.

QUESTION
I'm currently reviewing Mapbox's Android API for use as a backend tile provider, but the 6k tile limit is too small for our use. I saw this post regarding self-hosting Mapbox vector tiles and was wondering if that will allow me to bypass the tile cache ceiling limit.


texts: Getting back to side bar after detaching layer panel in QGIS

texts: I accidentally clicked on the detach icon on the layer side panel of QGIS and now it floats and I cannot get it back to the side bar where I want it.


texts: When trying to rasterize a SpatialLine .shp, I get a raster with points instead of lines?

texts:
EDIT:
I've solved the problem, the resolution was simply not high enough to properly show the lines in the plot. My code now looks like this and works fine:

r <- raster(ncol=1000, nrow=1000, extent(sa))
res(r) <- 20
rivers <- rasterize(lines, r)

I've been trying (for two days now), to rasterize a SpatialLines .shp file in R, but every time I load the new raster, I get points instead of lines. I've tried gdal_rasterize and different values and fields, but nothing seems to work. This is a croped version of what it should look like with my DEM underneath:

But instead, it looks like this:

This is my code:

sa <- raster("studyarea.tif")
lines <- readOGR("Gewässer_Uckermark/gewnet25_bb_a.shp")
r <- raster(ncol=1000, nrow=1000, extent(sa))
rivers <- rasterize(lines, r)
rivers[!(is.na(rivers))] <- 1

You can get the .shp file from here: https://geobroker.geobasis-bb.de/gbss.php?MODE=GetProductInformation&PRODUCTID=B9D461F1-99A1-4C10-97B4-9C36C0BD40B9


texts:
I have a csv table. Each record includes a start node (lat1,lon1) and an end node (lat2,lon2).
I have already imported the .csv into Postgresql 8.4. Now I would like set the geometry column. I tried to use this sql statement:

UPDATE mytable
SET the_geom = ST_GeomFromText('POINT(' || lat1 || ' ' || lon1 || ' ' ||
               lat2 || ' ' || lon2 || ')',4326);

but I receive this warning message:

WARNING:  OGC WKT expected, EWKT provided - use GeomFromEWKT() for this
            ERROR:  new row for relation "mytable" violates
            check constraint "enforce_dims_the_geom"

Do you have any idea on how to add geometries in my case?

Thank you
Selini


texts: Creating relationship database in ArcGIS Desktop

texts:
I am creating a file geodatabase with a feature class. I have created attribute columns in which a few of them have domain value YES/NO but when I select, NO, I have to provide a measurement value as per the condition.

How do I design such a database?

The objective of this database is to collect field information through the ArcGIS collector App where the field workers will measure the length of the roadside drainage if it does not satisfy the length specification.


texts: Create single map image from two images with different color tables

texts: Is it possible to merge the color tables of two images when using gdal to create a single image? Or is there another way to merge said images not using a color table?


texts: Calculating percentage of overlap between two different layers using ArcMap

texts: I have two different layers (Layer 1 and Layer 2) both made up of many thousands of polygons. Many of the polygons in Layer 1 overlap to differing degrees with Polygons in Layer 2. I need to calculate how much (percentage) each polygon in Layer 1 overlaps with the polygons in Layer 2.


texts:
Is there an easy way to convert MapInfo TAB files into GeoJSON?

I am working on a web-based mapping system which will use GeoJSON to transfer data. I also use MapInfo and would like to be able to convert files from MapInfo into GeoJSON.

I have tried the Universal Translator that comes with MapInfo but it doesn't support GeoJSON.

If there isn't a way to do it I will make a MapInfo plug-in myself (and probably post it online). I just don't want to waste my time making this if there is already one out there.

Update 2012-03-19: You can now export to GeoJSON through QGIS.


texts: PostGIS geometry to GEOS geometry without lwgeom dependency?

texts:
I've been looking for ways to convert PostGIS to GEOS geometries without linking to or dependency on lwgeom/liblwgeom.

The closest thing I can find is a mention of the pointcloud approach here:

  For an alternate example of an extension with a lighter dependency on
  PostGIS, check out pgpointcloud, which has it's own structure for
  spatial data (a point patch) and exchanges data with PostGIS via
  well-known-binary. This removes the liblwgeom dependency, which means
  it's possible to compile and use pgpointcloud without PostGIS
  installed, which is not entirely uncommon.

But I do not want to reinvent geometry types, and instead want to use what GEOS already have (, which I understand, corresponds to be subset of PostGIS/lwgeom types). So I am thinking about the path

PostGIS geometry  ==> (step1. via ST_AsEWKB or its like) EWKB
==> (step2. via GEOSGeomFromWKB_buf) GEOS geometry ==> call geos functions...

But I am not familiar with the code structures of the packages.

What I know so far are:

ST_AsEWKB for step1, but I guess I'll need a C/C++ equivalent.
Some examples reading WKB in step2.

/* WKB input */
g2 = GEOSGeomFromWKB_buf(uptr, size); free(uptr);

Can anyone with the expertise give some pointers to completing the loop?
(e.g. about specific pointcloud functions, better work flows, etc.)

(My target environment is a PostgreSQL extension (PostgreSQL 11) with GEOS 3.7 and PostGIS 2.5 as of Jan 2019.)

Other references:

How to use GEOS/C++ to efficiently find all point pairs closer than a threshold?


texts: OpenLayers 5: change the color of the highlighted feature

texts:
I've made a little script to highlight a feature on hover the mouse My inspiration comes from here. Now I would like to change the color of the highlighted feature. I can do this using style function with a vector layer but if I use ol.style.Style nothing change but the highlight run fine.

var highlightMouseOnHover = function() {
    selected_feature = new ol.interaction.Select({
      condition: ol.events.condition.pointerMove
    });
    map.addInteraction(selected_feature);
    selected_feature.on('select', function(e) {
      style: new ol.style.Style({
              stroke: new ol.style.Stroke({
                  color: 'rgba(121,121,125,1.0)',
                  lineDash: null,
                  lineCap: 'butt',
                  lineJoin: 'miter',
                  width: 0,
              }),
              fill: new ol.style.Fill({
                color: '#0d0887',
              }),
          });
    });
};

highlightMouseOnHover();

How I can solve this?
I'm not an expert of javascript


texts: PostGIS linestring to shapefile error with shapely/fiona

texts:
I have a PostGIS layer with linestring and I am trying to write the layer to shapefile using fiona and shapely.

Here is what I have so far:

from shapely.geometry import mapping,LineString
from shapely.wkt import loads,load
import fiona
from fiona.crs import from_epsg
import psycopg2
conn = psycopg2.connect("dbname='postgres',etc..) #connecting to DB
cur = conn.cursor()  #setting up connection cursor
crs = from_epsg(4269)
schema = {'geometry': 'LineString', 'properties': {'State': 'str', 'Date': 'str'}}
with fiona.open(r"routes.shp", "w", "ESRI Shapefile", schema, crs) as output:
    cur.execute("""select state,date,ST_AsText(geom) from trips;""")
    rows = cur.fetchall()
    conn.commit()
    for row in rows:
        print row
        g = loads(row[2])
        output.write({'properties':{'State': row[0], 'Date': row[1]},
        'geometry': mapping(g)})

It throws an error when it gets to:

('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)')
Traceback (most recent call last):
  File "C:\Users\Ralph\Documents\PythonScripts\LocationHistory.py", line 57, in <module>
    g = loads(row[2])
  File "C:\Python27\ArcGIS10.1\lib\site-packages\shapely\wkt.py", line 10, in loads
    return geom_from_wkt(data) #factory(geom)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\shapely\geometry\base.py", line 55, in geom_from_wkt
    "Could not create geometry because of errors while reading input."
ReadingError: Could not create geometry because of errors while reading input.

I used the same select query in postgres on the routes table and

('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)')

seems to be the error.

It writes the shapefile up to this record. what is strange is if I export the table to a shapefile using the PostGIS shapefile and dbf loader it exports the entire layer to a shapefile correctly.

My question is twofold:

How did postgis make a linestring from only one coordinate?
How would I go about circumventing this problem?

shapefile from PostGIS - this record ('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)')  which throws the error in shapely when trying to create a linestring (because it is only 1 coordinate) shows up in the PostGIS shapefile as a record with no geometric attribute


texts: Fix Processing AttributeError: 'NoneType' object has no attribute 'getCopy'

texts:
I have a problem with processing plugin of QGIS 2.13 and in particular with its Modeler.

When I have tried an operation, this is the log error:

*An error has occurred while executing Python code:
AttributeError: 'NoneType' object has no attribute 'getCopy'
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerDialog.py", line 424, in addAlgorithm
    self._addAlgorithm(alg.getCopy())
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerDialog.py", line 443, in _addAlgorithm
    self.repaintModel()
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerDialog.py", line 374, in repaintModel
    self.scene.paintModel(self.alg)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerScene.py", line 91, in paintModel
    item = ModelerGraphicItem(alg, model)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerGraphicItem.py", line 60, in __init__
    self.pixmap = element.algorithm.getIcon().pixmap(15, 15)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\modeler\ModelerAlgorithm.py", line 121, in algorithm
    self._algInstance = algList.getAlgorithm(self.consoleName).getCopy()
AttributeError: 'NoneType' object has no attribute 'getCopy'
Versione Python: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
Versione di QGIS: 2.18.13 Las Palmas, 8cba0bb9eb
***Path di Python:
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python
C:/Users/Giuseppe Cillis/.qgis2/python
C:/Users/Giuseppe Cillis/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.18\bin
C:\PROGRA~1\QGIS2~1.18\apps\Python27
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/Giuseppe Cillis/.qgis2//python*

I tried to unistall and install and I also tried to delete .qgis2 folder .

How can I solve this problem?


texts: GeoServer NetCDF loader ignores the scaling factor - how to fix?

texts:
I am using GeoServer 2.13.2 on an Ubuntu Trusty 14.04 Linux docker container and I have most plugins installed and certainly the NetCDF plugin.

For the most part it seems work well and I have all layers exposed and working until I discovered that the values were completely wrong when queried.

I just used the Layer Preview to test the values. The values were gigantic. After running toolsUI-4.6.jar and examining the structure I see:

netcdf
 /data/weather/June2018_WindWaveSwellSSTSSPUV1.nc {
  dimensions:
    longitude = 480;
    latitude = 241;
    time = UNLIMITED;   // (241 currently)
  variables:
    float longitude(longitude=480);
      :units = "degrees_east";
      :long_name = "longitude";

    float latitude(latitude=241);
      :units = "degrees_north";
      :long_name = "latitude";

    int time(time=241);
      :units = "hours since 1900-01-01 00:00:0.0";
      :long_name = "time";
      :calendar = "gregorian";

    short sst(time=241, latitude=241, longitude=480);
      :**scale_factor = 5.742090041944516E-4**; // double
      :add_offset = 287.4019345751854; // double
      :_FillValue = -32767S; // short
      :missing_value = -32767S; // short
      :units = "K";
      :long_name = "Sea surface temperature";

    short sp(time=241, latitude=241, longitude=480);
      :**scale_factor = 0.788649718939313**; // double
      :add_offset = 78192.71700326553; // double
      :_FillValue = -32767S; // short
      :missing_value = -32767S; // short
      :units = "Pa";
      :long_name = "Surface pressure";
      :standard_name = "surface_air_pressure";

    short u10(time=241, latitude=241, longitude=480);
      :scale_factor = 0.001169629851026868; // double
      :add_offset = 3.466105248551049; // double
      :_FillValue = -32767S; // short
      :missing_value = -32767S; // short
      :units = "m s**-1";
      :long_name = "10 metre U wind component";

It would seem the scaling factor is being ignored.
Is there a way to add it to the index schema file or how can I correct it so the layer returns the correct values?


texts: Determine shortest distance between point to a (small) path of OSM?

texts:
When using offline on a mobile (Android) device: How can I calculate the shortest distances between any point and the nearest (small) OSM path?

The offline file with graph/routing data could be used by a library to do the distance calculating. The size could be for a country like the Netherlands.

What kind of paths? I mean especially the smallest paths, like footpaths, hiking trails, etc. I don't need any further routing facilities.

How will I use the data? About 6 years I do geocaching with my onw geocaching app. I use OsmDroid, OsmDroid-with-Mapsforge i.c.m. OsmBonuspack. I work with downloaded mapsforge-v5 downloaded maps.

I am building a lot of AI stuff in my app so that I can guesstimate locations when I can't figure them out. For the last piece of AI I need to determine the shortest distance between any point to the nearest (foot)path.

If offline is too difficult, is there any API that can be called?

An example: I have a series coordinates (or points) in the woods. Which point is nearest to a footpath?

UPDATE: a lot of useful info is found at this location. So, now I only need to use the unpacked and pre-processed stuff to use the detailed API. Anyone experience with a simple main() method that does the trick?


texts: Simulation in grid/array over map in Python

texts:
I'm new to GIS.
I'm a physicist and I have a simulation which I usually run in numpy arrays (usually ending up with a heatmap or contour plot like below).

However, I now want to run the simulation over the geography of a country or region rather than just a square. I have the shape file for the region I'm interested in, and I know how to plot basic lat/long points, but how can I create some sort of array or grid to run a simulation over?
I've tried turning the shapefile into a binary mask and putting that into a numpy array, but that gets rid of any spatial information.

I'm solving a differential equation over a discrete scalar field (i.e. my numpy array is a discrete grid over some region, and the values in the array are evolved over time). I'm interested in how the geometry of the system affects the answer. Originally I used a square array and experimented with simple binary masks like circles etc., but the model is broadly to do with population movement so I want to use the geometry of a real region to investigate things further.


texts: Importing large amounts of CADD data into ArcMap

texts:
I work at a utility coordination company and although I have a background in coding and computers I'm brand new to ArcMap and GIS so I apologize if I come across as a little oblivious.

Basically, I am trying to import lot of of CADD data (around 2 GB) into ArcMap and hopefully turn it into a geodatabase or shapefile. But, when I try to import all this data into ArcMap, it just freezes up and lags very badly.

This is an issue I have been having in general with ArcMap and I suspect updating the company computer to Windows 10 will update some of the drivers and potentially help with this, but I am curious if there is a better way?

This is a lot of data, but I am sure this has been done before successfully.


texts: Exporting DXF contour lines in QGIS. Units Problem

texts:
I am trying to export the contour lines of a topography from QGIS to DXF AutoCAD to generate a 3d model of it. The data comes from a source that has the elevation in FT. It matches with my coordinates system in plan but when I export it the Z Height remains in ft. How can I change that in QGIS so I get my model in meters?

Still having the same issue. Once I am in Rhino or Cad, elevation data is at scale but then the distances horizontally do not match at all (are scaled more than 10.000 times).
Sorry to be such a pill... The file I am using to be more specific is this one https://geohub.lacity.org/datasets/contour-lines-2006.


texts: Different buffer distance left and right form line?

texts: I'm trying to make different buffer zones, left and right from the line. I made processing model (picture). I can do this with this module, but I can not specify the value that I want. I am executing only the default value in the model. Is there a possibility that, when starting the module, open the window to enter the value,  left and right?


texts: Calling OpenWeatherMap API in Leaflet.js map

texts:
How to get weather data from OpenWeatherMap API in a popup by clicking any location on the map?

i am still new at calling APIs in my map, API return JSON files but i don't know how to call the API in my map and how to bind JSON file in an HTML element.

This is the API which uses coordinates to return weather data in a JSON file: api.openweathermap.org/data/2.5/weather?lat=35&lon=139

but i need to replace these lat and lng parameters by clicking any location on the map.


texts: Find the area of a map (Xkm x Ykm) from the left & right longitude and latitude

texts:
How can the area of a map from OpenStreetMap be determined in the form of Xkm x Ykm, given the following information:

Left Longitude  = 103.554879
Right Longitude = 103.740497
Top Latitude    = 1.585770
Bottom Latitude = 1.490873

Any ideas?


texts: Extract: error writing to file qgis-bin.exe

texts: I have been facing problem during the installation process. The installation is getting aborted by showing this error message. Since this happened with 2.18 version, I tried with 2.16.3 version. Even then i am getting the same message. Please help me in overcoming this problem.


texts: .tif to PCraster format conversion in QGIS 3.2

texts:
I am using QGIS 3.2 and trying to convert a resampled .tif file into PCRaster format which is used for a downstream model.

Unfortunately when i try to export the resampled.tif file into PCRaster format i get the following error:
"Unsupported Data Source: resampled.map is not a supported raster data source"

Can anyone help to export my tif to PCRaster?


texts: Combining two TAB files in MapBasic - losing data

texts:
I am using mapbasic to perform some calculations on a table. Everything is fine until I combine two .tab files.

Both files have identical table structures, however upon combination some numerical (float) values are being set as 0.

Values from HUNDRED_CL:

Values after combining the files:

The CR1 data is fine, but as soon as I append the CL1 data to it, some TO values become 0, can somebody assist?

Table structures for reference:


texts: How to plot latitude, longitude and elevaltion data in python

texts:
I have latitude and longitude and elevation data in a CSV file and I want to process it in 3d to something like this

the data set is kind of like this

I just want latitude, longitude and elevation


texts:
I have some code in a function that uses a list, which works, but is there a way to write a list outside of the function and use it within the function code?

For example in my code below I have a list called variables and then each type in variables is again listed in the WOTable with the code remap['type'] but I would like a code that does not relist the variable again with remap[type] in WOTable rather have a list outside the function where the list parameters can change but the function code stay the same. This script is for a weighted overlay in ArcGIS 10.2.

I have tried calling the list as weighted_overlay(variables, other parameters), and using just variables = [] in the function code, which did not work.

def weighted_overlay(raster1, raster2, raster3, out_raster):
    in_raster1 = arcpy.Raster(raster1)
    in_raster2 = arcpy.Raster(raster2)
    in_raster3 = arcpy.Raster(raster3)

    variables = ['remap1', 'remap2', 'remap3']
    remap = {}
    for type in variables:
        remap[type] = RemapValue([[1, 1], [2, 2], [3, 3], [4,4], ["NODATA", "NODATA"]])

    WOTable = WOTable([[in_raster1, 50, "VALUE", remap['remap1']],
                   [in_raster2, 35, "VALUE", remap['remap2']],
                   [in_raster3, 15, "VALUE", remap['remap3']]], [1, 9, 1])

    outWeightedOverlay = WeightedOverlay(WOTable)
    outWeightedOverlay.save(out_raster)
    return

Sorry if this has been asked before but I was unable to find a solution and I do not have a lot of experience with python lists. Eventually I would like to make this a class function as part of a larger analysis script to reuse.


texts: SQL carto.js API getting the sum of a variable

texts:
I need to sum up the values of a given field rather than counting . The code creates the variable 'sum' that gives you the number of rows, what if need the sum of the values of a given field/column?? is there any sum() operator I can include to get it?

//sql API to filter and then js to count the number of rows filtered

sql.execute("SELECT * FROM table_4mi_dataset_p5 where type in ("+ incident + ") and origin in ("+ origin + ") and sex in ("+ gender + ")")
.done(function(data){
 for(i = 0; i < data.total_rows; i++){ // loop that reads each column of the table
    sum = data.total_rows}

console.log(sum)

});


texts: Python shapely difference between two polygons

texts:
I have a geometry question using shapely,Fiona and polygons shapefiles with big data(over 700000 features any polygon shapefile).

polygon1.difference(polygon2)

I want to find the difference between two polygons using shapely (or maybe geopandas?) and I have read this post, I want the difference calculate only in the polygons segments(lines) and not the polygon geometry.

for example if I use code from the post with add more lines the final export difference is a polygon layer,I want the line layer to show me the differences between two polygons.

with the light blue line is the difference.

that because I think is more fast the difference calculator to create lines from polygons.


texts:
I am trying to use Voronoi polygons algorithm in QGIS in a dataset of points but strangely I am getting an overlap of polygons.
These points correspond to the vertices of a set of polygons.

When using MMQGIS the polygons are correctly defined, however I can't use this plugin for the number of points I need to use (much higher number of points) and since I can't set a buffer in this algorithm.

Attached you will find some screenshots of the results I got:
- Fig1 is the set of polygons from which were extracted the vertices
- Fig2 are the voronoi polygons obtained from the native qgis algorithm
- Fig3 are the voronoi polygons obtained from MMQGIS plugin

It is observed from the selection that there is an overlap of some polygons.

How can I avoid this?

Figure 1 - Polygons from which were extracted the vertices

Figure 2 - Voronoi output from QGIS native algorithm

Figure 3 - Voronoi output using MMQGIS plugin

Dataset of Points


texts: ESRI JavaScript API FeatureLayer Click Event

texts:
I am using FeatureLayers in the API to load in a layer.  I have an info template set for the layer for identifying individual features.

I am trying to disconnect the click event from the layer while I am drawing a polygon for another tool.  However, I haven't found a good way to disconnect that handler since it isn't one I assign.

I have tried layer.disableMapEvents() which didn't work.  I don't have an info template set for the map itself either, so I don't think that one is firing.

I have searched through some other posts here, and so far the only solution I have found is to attach another handler to the map click event and call map.infoWindow.hide() while I am doing the drawing.  This still goes through the identify process but keeps the window away.

I also tried calling the map.graphics.disableMouseEvents(), but that doesn't work either.

I am looking for a way to either handle the map click event myself and disconnect ESRI's standard way for the FeatureLayer.  I just want to be able to shut off identifying of features temporarily.

A way to disconnect the event and hook it up again would be great.  That would help if I want to attach it to another mouse event, too.


texts: Flow modeling and FOSS GIS; lahars and HEC-RAS

texts:
What are the options for FOSS flow modeling of dilute debris flows in QGIS or GRASS?
Have you, or can you refer me to studies where the authors, tested these models/plugins/functions/etc. for use in lahar modeling?
  (Lahars are muddy volcanic flows, I'm interested in ones that will have a very small solids fraction of <10-20%, making some flood models applicable to the study; other add'l info: my field site has poor DEM coverage and is highly channelized, aka not a flood plain)

I am working on developing a FOSS workflow for disaster risk reduction, specifically for using DEMs in FOSS GIS that can either be used as input in separate flow models or with flow models that run in GRASS or QGIS. I know of a few functions that can do this (r.inund.fluv, r.hazard.flood, r.drain, etc) but they seem to be untested for this function.

I was interested in using HEC-RAS but my collaborator feels that it only runs effectively in Arc as HEC-geoRAS. The other model I was interested in (again, Arc) is the D8-algo based LAHARZ model, which (it seems?) is kind of similar to r.drain, at least in theory.

Any information at all on the matter is greatly appreciated!


texts:
I am trying to fill a field and I am using arcpy's CalculateField_management function.

I want to fill the new field with the name of the shapefile, so I am using arcpy.Describe. But when I run the script I get this error: ExecuteError: ERROR 000539: Error running expression: PRUEBAS Traceback (most recent call last): File "", line 1, in  NameError: name 'PRUEBAS' is not defined. 'Pruebas' is the name of the shapefile, I obtain it from desc.baseName

The code is:

import arcpy, os, sys
from arcpy import env
from arcpy import Describe

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "FIELD"
        self.description = ""
        self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""

# First parameter

    param0 = arcpy.Parameter(
        displayName="Shape entrada",
        name="in_features",
        datatype="DEShapefile",
            parameterType="Required",
        direction="Input")

# Second parameter

param1 = arcpy.Parameter(
        displayName="Nuevo campo",
            name="field",
        datatype="Field",
            parameterType="Required",
        direction="Input")

    parameters = [param0, param1]
    return parameters

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True
def updateParameters(self, parameters):return
def updateMessages(self, parameters):return
def execute(self, parameters, messages):
    """The source code of the tool."""
    a= parameters[0].valueAsText
    b= parameters[1].valueAsText

    desc = arcpy.Describe(parameters[0])
    a=arcpy.AddField_management(a, b, "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

    a=arcpy.CalculateField_management(a, b, desc.baseName,"PYTHON_9.3",code_block="")

    return a

`


texts:
I have animal relocation data and want to know how many times that animal walked within a polygon for more than one step (i.e there is more than one point within a polygon). I am more interested in how many times this happens, and not necessarily how many points there were consecutively (I suspect that there will only be 2 relocations in a row in the polygon, or maybe 3).

I don't want to count the number of points in the polygon because it needs to be: In the polygon, in the polygon... and NOT in the polygon, out of the polygon, in the polygon.

Any ideas?

This question still needs to be answered! Help please! I've attached an image that shows when I do want to count the path (red), and when I don't (green).

Note: I want to count how many times this (red in the picture) happens, not how many points there are.


texts: Digitized GIS Data - Wrong Coordinate System, Right Location

texts:
I (correctly) georeferenced a JPEG to an existing shapefile that was in the correct coordinate system.

I then proceeded to create line feature classes by manually tracing some of the features (roads) from the JPEG.

The mistake: When creating the line feature classes, I used the wrong coordinate system for the feature datasets, and so the line feature classes are all in the wrong coordinate system.

The result is that the feature classes are appearing in the correct location (because I manually traced them over the georeferenced JPG) but are in the incorrect coordinate system. Essentially I need a way to update the coordinate system to the correct projection without changing the location of the lines.


texts: Changing displayed elements of a KML file based on current system time in Google Earth

texts:
Is there any way to have Google Earth turn features in a KML file on and off automatically based on the current system time? I have created the simple example below of a rectangle with a <TimeStamp> set to end before now, but it displays anyway when I load the file. If I set it for a minute or so in the future then load it and wait, the rectangle doesn't disappear when the system clock rolls past the specified end time.

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<kml xmlns="http://earth.google.com/kml/2.2">
<Document>
  <Placemark>
    <TimeSpan>
      <end>2013-07-30T22:38:00Z</end>
    </TimeSpan>
    <Polygon>
      <extrude>0</extrude>
      <tessellate>1</tessellate>
      <altitudeMode>clampToGround</altitudeMode>
      <outerBoundaryIs>
        <LinearRing>
          <coordinates>
            -10,10
            -10,-10
            10,-10
            10,10
            -10,10
          </coordinates>
        </LinearRing>
      </outerBoundaryIs>
    </Polygon>
  </Placemark>
</Document>
</kml>


texts: Convert a List to a file list to use in an arcpy tool

texts:
I just can't figure out how to feed a list of files into an arcpy tool...

I have

 List_tif.append(os.getcwd()+"\\Stage1\\"+name+"\\GeoTif\\"+name+"_"+str(ddp.pageRow.PageName) + ".tif;")

That generates the list and then I try to feed it into arcpy

arcpy.MosaicToNewRaster_management(List_tif,os.getcwd()+"\\Stage2\\jp2", name+"_tif.jp2", "90019.prj", "", "", "3", "LAST","FIRST")

but I get an error as the tools is expecting file1.tif; file2.tif etc but is getting a List

  'P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\Stage1\LS\GeoTif\LS_O10.png;',
  'P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\Stage1\LS\GeoTif\LS_O11.png;',
  'P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\Stage1\LS\GeoTif\LS_O12.png;',
  'P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\Stage1\LS\GeoTif\LS_O13.png;']

I tried

PNGList=string.split(List_png, ";")

but this doesn't work either. This should be simple (:-

EDIT 1:

Thanks guys...I tried the list rasters as well but it doesn't even generate the list

>>> rasters=arcpy.ListRasters('*')
>>> rasters
>>> files=arcpy.ListFiles()
>>> files
>>> import glob
>>> rasters=glob.glob('*.png')
>>> rasters
['BA25k_A2.png', 'BA25k_A3.png', 'BA25k_A4.png', 'BA25k_A5.png', 'BA25k_A6.png', 'BA25k_A7.png', 'BA25k_A8.png', 'BA25k_B10.png', 'BA25k_B11.png', 'BA25k_B2.png', 'BA25k_B3.png', 'BA25k_B4.png', 'BA25k_B5.png', 'BA25k_B6.png', 'BA25k_B7.png', 'BA25k_B8.png', 'BA25k_B9.png', 'BA25k_C10.png', 'BA25k_C11.png', 'BA25k_C12.png', 'BA25k_C2.png', 'BA25k_C3.png', 'BA25k_C4.png', 'BA25k_C5.png', 'BA25k_C6.png', 'BA25k_C7.png', 'BA25k_C8.png', 'BA25k_C9.png', 'BA25k_D10.png', 'BA25k_D11.png', 'BA25k_D12.png', 'BA25k_D13.png', 'BA25k_D2.png', 'BA25k_D3.png', 'BA25k_D4.png', 'BA25k_D5.png', 'BA25k_D6.png']
>>>

I didn't mean to ask this specifically for the arcpy mosaic tool as passing the list of files is important for lots of other tools as well such as creating a mosaic dataset, filling a gdb etc.

I also tried your vtab example

vtab = arcpy.ValueTable(1)
for tif in List_tif:
    vtab.addRow(tif)

print "Generating jp2 from png tiles @"+str(datetime.datetime.now())
##PNGList=string.split(List_png, ";")
arcpy.MosaicToNewRaster_management(vtab,os.getcwd()+"\\Stage2\\jp2", name+"_tif.jp2", "90019.prj", "", "", "3", "LAST","FIRST")

but I still get the bands error.


texts: Running Python IDLE only works with remote into desktop?

texts:
An extremely odd issue.  I have python scripts I use, pretty straight forward.  I use them to automate printing of maps in Arcmap (10.2).  I work from home and in my office.  When I work from home I remote into my desktop at my office. I open IDLE to run one of my scripts, If I remote into my computer it works fine.  When I am in the office and using my desktop directly, IDLE crashes saying "Python.exe has stopped working" and that's it.

Any ideas what could be causing this or potential fixes?

I am running python 2.7


texts: Converting MapBox PBF vector tiles to GeoJSON?

texts:
There are many pieces MapBox in file.vector.pbf. I can get information from them in GeoJSON using tippecanoe, but the geometry will be broken(screen below).

How can I normalize geometry or combine all the pieces into one big card?


texts:
I have calculated atmospherically corrected reflectance values for Landsat 8.

From that, I am getting correct values of NDVI in QGIS, but for EVI I am not getting the correct values. For some images it goes beyond 1 and for others it shows black output in QGIS.

I am using the EVI formula derived for MODIS.


texts: Data Converting into Collection Tileset Instead of Type (Line) Tileset in Mapbox Studio

texts:
A recent issue I've encountered with Mapbox Studio is the apparent inability to upload a geojson as anything other than a collection tileset, instead of a point, line, or polygon tileset. For example, the old tileset (circled in black) shows the line icon where the new one (circled in blue) shows the collection icon.

My difficulty lies in updating the style. When I go to switch out the data Mapbox Studio detects the different data type and resets the style...which means I have to style that layer by hand each time I update the information. This becomes time consuming when the tileset is used on 15+ layers. So far Mapbox's response has been to "update the style" which is not helpful as that's what I'm trying to do.

I put the geojson into a text editor to remove "feature collection" and just make it a list of features, but Mapbox won't upload it and errors out. I created a dataset directly in Studio but when exported to a tileset the same thing happens (becomes a collection instead of line). I have also tried this with data that contained a single polygon or a single point and encounter the same issue.

Does anyone have any suggestions? I tried adjusting the data so maybe the issue is in learning how to save the styling information or preventing the style from clearing when a different data type is detected?


texts: I want to create a buffer in meters around points in Europe. But I cannot find a proper CRS for Europe that is in meters instead of degrees?


texts: Can't get the correct version of GDAL installed to use ReadAsArray without getting error in Python

texts:
Using the instructions here . . . https://pythongisandstuff.wordpress.com/2011/07/07/installing-gdal-and-ogr-for-python-on-windows/  I downloaded the following files:
GDAL-1.11.1.win-amd64-py2.7.msi and
gdal-111-1800-x64-core.msi

After doing this I was able to create my Python script (which manipulates DTED data), and run this script from my C# code.

Then I ran into a few problems when attempting to manipulate the pixel values of the DTED file:
1)  I'm trying to use ReadAsArray and I get "ImportError: No module named _gdal_array" on the line in my Python script . .  data = rb.ReadAsArray(0, 0, cols, rows)
Referred to:
Cannot import gdal_array
http://lists.osgeo.org/pipermail/gdal-dev/2011-May/028708.html
2)  I also thought that I could use gdallocationinfo . .  http://www.gdal.org/gdallocationinfo.html
I get a syntax error on the file name . . .

      gdallocationinfo C:\Data\TestFiles\TST.DT1 1000 1000
      SyntaxError: invalid syntax

      I find that it gives an error message on any file name that I enter.

Therefore, I installed OSGeo4W in order to be able to use some of the above mentioned GDAL functions that I would need (after reading this https://hub.qgis.org/issues/8811).
By running the OSGeo4W Shell, I was able to type this in at the command line and get the desired result:

C:>gdallocationinfo C:\Data\TestFiles\TST.DT1
 1000 1000
Report:
  Location: (1000P,1000L)
  Band 1:
    Value: 859

The problem is that I don't think that I can use a Python script with this OSGeo4W shell, which means that I would have to run some other kind of script from my C# code, I guess.

Are there other versions of msi files that I should have downloaded instead of GDAL-1.11.1.win-amd64-py2.7.msi and gdal-111-1800-x64-core.msi?

I'm just not sure which approach I should take.


texts:
I can draw a polygon in OpenLayers by running this code.

        var points = [];

        var point = new OpenLayers.Geometry.Point(0, 0);
        point.transform('EPSG:4326', 'EPSG:900913');
        points.push(point);

        var point = new OpenLayers.Geometry.Point(0, 30);
        point.transform('EPSG:4326', 'EPSG:900913');
        points.push(point);

        var point = new OpenLayers.Geometry.Point(30, 30);
        point.transform('EPSG:4326', 'EPSG:900913');
        points.push(point);

        var point = new OpenLayers.Geometry.Point(30, 0);
        point.transform('EPSG:4326', 'EPSG:900913');
        points.push(point);

        var ring = new OpenLayers.Geometry.LinearRing(points);
        var polygon = new OpenLayers.Geometry.Polygon([ring]);

        // create some attributes for the feature
        var attributes = { name: "my name", bar: "foo" };

        var feature = new OpenLayers.Feature.Vector(polygon, attributes);
        select.addFeatures([feature]);

        map.addLayer(select);

My map can show this polygon perfectly. Here I did point.transform to match my map's projection system.

Now I want to edit this polygon like 'DragFeature', 'SelectFeature', 'Separator', 'DrawHole', 'ModifyFeature' etc by OpenLayers.

How can I do that?


texts: Difficulty joining layers in ARCGIS: Always asked to define projection

texts:
Having followed manual instructions I am experiencing problems joining data to a map. What I have is a UK Districts shapefile and imported data with my coordinate data in it.

This coordinate data I produced in excel by assigning centroids to each district.
On advice I imported my excel file as a Personal geodatabase(.mdb) instead because apparently(ArcGis Join Data CSV Problem) ARCGIS treats this better but I still cannot fully merge layers.

The coordinates to 'map' onto my map  correctly - see image - but the two layers are not joined. I have the error that I need to define a projection. I try to do this through the ARC Toolbox but it still says unknown coordinate system in properties. I then add Join/Relate from the properties window but get another error message.

What Can I do to diagnose the problem here? I've been stuck for so long!


texts:
I'm trying to use GDAL in python to conditionally select raster bands from a grib file based on the time when the weather feature was measured (think temperature).

If I run GetMetadata() on the dataset as a whole, I get no information.

path = r'...\ds.snow.bin'
dataset = gdal.Open(path)
dataset.GetMetadata()

  {}

However, I do get information if I run GetMetadata() on a specific band:

dataset.GetRasterBand(1).GetMetadata()

  {'GRIB_COMMENT': 'Total snowfall [m]',  'GRIB_DISCIPLINE':
  '0(Meteorological)',  'GRIB_ELEMENT': 'SnowAmt',
  'GRIB_FORECAST_SECONDS': '0 sec',  'GRIB_IDS': 'CENTER=8(US-NWSTG)
  MASTER_TABLE=1 LOCAL_TABLE=0 SIGNF_REF_TIME=1(Start_of_Forecast)
  REF_TIME=2019-05-29T15:00:00Z PROD_STATUS=1(Operational_test)
  TYPE=1(Forecast)',  'GRIB_PDS_PDTN': '8',
  'GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES': '1 29 2 0 0 255 255 1 0 1 0 0
  255 -1 -2147483647 2019 5 29 18 0 0 1 0 1 2 1 6 1 0',
  'GRIB_PDS_TEMPLATE_NUMBERS': '1 29 2 0 0 0 255 255 1 0 0 0 0 1 0 0 0 0
  0 255 129 255 255 255 255 7 227 5 29 18 0 0 1 0 0 0 0 1 2 1 0 0 0 6 1
  0 0 0 0',  'GRIB_REF_TIME': '  1559142000 sec UTC',
  'GRIB_SHORT_NAME': '0-SFC',  'GRIB_UNIT': '[m]',  'GRIB_VALID_TIME': '
  1559152800 sec UTC'}

Is there a way to get this metadata for the entire grib file at once so I can parse over it and select only bands that fall within a specified time?

I know cfgrib using XAarray can do this, but I would rather use GDAL.


texts: Geoserver 2.2.5 tomcat unable to load servlet2.5.jar

texts:
am trying to use geoserver 2.2.5 with tomcat 7.0.34 following this tuto:
http://ian01.geog.psu.edu/geoserver_docs/software/geoserver.html
but am getting this error concerning the servlet jar:

     5 mai 2013 14:10:25 org.apache.catalina.core.StandardEngine startInternal
     INFO: Starting Servlet Engine: Apache Tomcat/7.0.34
     5 mai 2013 14:10:25 org.apache.catalina.startup.HostConfig deployWAR
     INFO: Déploiement de l'archive C:\Users\OUK\Desktop\PFE\tomcat-7.0.34\apache-tomcat-7.0.34\webapps\geoserver.war de l'application web
     5 mai 2013 14:10:26 org.apache.catalina.loader.WebappClassLoader validateJarFile
     INFO: validateJarFile(C:\Users\OUK\Desktop\PFE\tomcat-7.0.34\apache-tomcat-7.0.34\webapps\geoserver\WEB-INF\lib\servlet-api-2.5.jar) - jar not loaded. See Servlet Spec 2.3, section 9.7.2. Offending class: javax/servlet/Servlet.class
     log4j:WARN File option not set for appender [geoserverlogfile].
     log4j:WARN Are you using FileAppender instead of ConsoleAppender?
     05 mai 14:10:35 WARN [config.CustomEditorConfigurer] - Passing PropertyEditor instances into CustomEditorConfigurer is deprecated: use PropertyEditorRegistrars or PropertyEditor class names instead. Offending key [org.geotools.util.Version; offending editor instance: org.geoserver.platform.util.VersionPropertyEditor@126456a
     5 mai 2013 14:10:35 org.geoserver.platform.GeoServerExtensions checkContext
     GRAVE: Extension lookup occured, but ApplicationContext is unset.
     05 mai 14:10:35 ERROR [geoserver.global] -

can anyone figure out how to solve this problem, thnks.


texts: Nearest distance from point to pixel in QGIS?

texts:
I have 1) a point layer and 2) a binary raster file (1 and 0). I need to estimate the Euclidean distance from each point to the nearest "1" pixel. I know there is a tool in ArcGIS but I only have access to QGIS/GRASS; in GRASS there is a clear way to estimate this distance but it's between two vectorial layers (using v.distance).

So far, I tried:

1) in QGIS, the Proximity tool, but it doesn't give me the result I want, since it considers only one raster layer

2) in GRASS, r.distance but I get an error ("Raster map  is not CELL"), and I am still trying to figure it out, but I'm still not sure which values will it give me since I don't see any option to select a set of pixels ("1", in this case), instead of all the layer

Does anyone know how to do that in QGIS?


texts: How to make MapInfo able to edit geographic data from PostGIS?

texts:
I'm using MapInfo 10.5 to read geographic data from PostGIS 1.5.2. After a pg_dump/pg_restore to another PostGIS-database on another server I'm trying the same via a VPN-tunnel. When I open a table I get the message "Map objects cannot be edited" and then I can't edit the geometries. The other attribute data is editable.

I've tried to create my database and tables both by hand and with easyloader with the same result. At the same time when I restore the database on my local computer I only get the message for one of the tables.

Any idea what could causing this problem?

Regards Ulrik


texts: Stripes in Landsat-image after topographic correction (SAGA)

texts:
I got stripes in my topographic corrected Landsat images (see second image below).
I used SAGA GIS 2.0.8 at xubuntu 12.04

Image and DEM are in the same projection, tried it with  EPSG:32636 and EPSG:4326


texts: Save shapefile as SQL script with ogr2ogr?

texts:
I have a some shapefiles and I want to convert them into SQL scripts.
I know that ogr2ogr can put data into a database itself. But I want to have a SQL script which can be executed and put data to DB.

Any advice?

If it impossible with ogr2ogr maybe I can use other tools for this?


texts: Intersection Polygon-PolyLine with pyqgis

texts:
I want to make intersection between polygon and polyline shapefile with pyqgis, but I get an empty output. Here the code that i use :

params = {'INPUT': polygone.shp, 'OVERLAY': line.shp, 'INPUT_FIELDS': [], 'OVERLAY_FIELDS': [], 'OUTPUT': 'memory'}
processing.run('native:intersection', params)

bellow the picture of polygon and polyline shapefile. The both shapefiles have the same crs : epsg:2154


texts:
Can I edit this code below to show only point layers?
This is being used in a Plugin for QGIS 2.
I am creating if that makes any difference.
self.addParameter(ParameterVector(self.INPUT, self.tr('Input layer'), [ParameterVector.VECTOR_TYPE_ANY]))


texts: How to change the pixel value of a raster dataset which overlaps a polygon featured shapefile in ArcGIS

texts: I have a raster data of 1992 and a 'shapefile' of Los Angeles county. I want to reclassify the pixel value of the raster data set to 'no data' or zero in the region where the raster and the 'shapefile' overlaps. I'm using ArcGIS 10.1. I tried using spatial analyst tool> 'reclass'. But I'm not getting the desired result.


texts: Increase spatial query performance (st_difference)

texts:
I have two link tables

one is 9029 rows (and let's say this table is A) and the other is 11698 rows (and this table is B)

the both tables are road link and A is for tunnel link and B is road link

From table road link (table B), I need to erase segment where tunnels are (Table A)

Basically, it is Table B - Table A.

This is my spatial query though, it took me so much time and I do not know why.. It seems like other GIS specialists also experiences that st_difference takes quite lots of time.

CREATE TABLE rd_tunnel_erased AS
SELECT b.objectid, b.roadclass, b.linkclass, b.linktype, st_difference (b.shape, a.shape) as shape
FROM b, a

As I said b is road link and a is tunnel link (and the two table does have same schema but difference geometry)

Here are two screenshots that might help you

this is overlayed a and b table

from b table, I erased (differenced) a table.

With only 1 link, the st_difference worked fine and fast.

Somehow, when it goes over my real data, it takes more than 3~4 days (in fact, I do not know how many days and hours it took.. because it failed for trivial mistakes. I am afraid that I might have to waste my time for re-running this query and ended up getting nothing.)

All my tables are spatially indexed.


texts: Using removeall map layers for QgsProject in PyQGIS?

texts:
I want to rewrite a projekt with new layers while owerwriting all old layers. I found the function removeAllMapLayers in the qgis.core api

But when i try:

project = QgsProject.instance()
project.removeAllMapLayers()

Where I am in a project with some layers. After trying the code I get this error:

Traceback (most recent call last):
File "<input>", line 1, in <module>
AttributeError: 'QgsProject' object has no attribute 'removeAllMapLayers'

Is there any other way to empty all layer in a project? I want to write to the same project when it's empty.


texts: Reading in GeoTIFF with different CRS and Transform with rasterio

texts:
I have multiple TIF files with different resolutions that I want to combine using rasterio. I am using a modification of the example given in the documentation here. However, on closer inspection whilst all the TIFs have the same CRS, they don't have the same affine transformations.
OPTION 1 Here is the code I have now:
    import rasterio as rio
    import numpy as np
    import numpy.ma as ma

    def tif_load(tif_path, dimensions, resampling_method='bilinear'):
        with rio.open(tif_path) as raster:
        # resample data to target shape
        data = raster.read(1, masked=True, out_shape=(dimensions[0], dimensions[1]),
                           resampling=Resampling[resampling_method]
                           )
        meta = raster.profile
        # scale image transform
        transform = raster.transform * raster.transform.scale(
            (raster.width / data.shape[-1]),
            (raster.height / data.shape[-2])
        )
       return data, meta, transform

I guess I could use the code above, but is there anyway to read in the tif file with a particular crs and transform? I was thinking I could use one of my TIFs to set a standard crs and transform and then use the above, passing in this standard crs and transform to read in the others? I would probably have to be careful to remember to modify the meta for the raster to make sure it reflects the destination crs and transform.
OPTION 2 Alternatively, I could read in the tif files like so
import rasterio as rio
import numpy as np
import numpy.ma as ma

with rio.open(tif_path) as raster:
   # resample data to target shape
   data = raster.read(1, masked=True)
   meta = raster.profile

And then I could set a standard crs and transform from one of my tifs and use reproject like they do in the documentation here. For example, say I set the tif_path as our standard crs and transform and I read in tif_path2 (source) with data2 and meta2, then we can do the following
src_transform=meta2['transform']
src_crs=meta2['crs']
source=data2
destination=ma.copy(data)
dst_transform=meta['transform']
dst_crs=meta['crs']
with rio.Env():
    reproject(
    source,
    destination,
    src_transform=src_transform,
    src_crs=src_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    resampling=Resampling.bilinear)

# Assert that the destination is only partly filled.
assert destination.any()
assert not destination.all()

Will these give the same result? Which is better? If this changes things I am reading in the data as a masked array.


texts: How to migrate PostGIS databases between Postgres versions?

texts:
Mac OS

I am trying to upgrade my Postgres version 9.2 to 9.6 and failing.  My strategy was to follow the standard advice for Postgres databases and use pg_dumpall.  This does not appear to work if you are using PostGIS extensions because the dump contains 1000s of references to the actual PostGIS library on disk.  See this.

So how does one move around PostGIS databases without links to actual libraries embedded in them?


texts:
In a QGIS map legend, what methods can be used to show examples of a feature's label, in addition to its symbol?

For example, a single point symbol might be used to show all cities, but the label style may be varied in size, color, typeface (bold, italic, underline etc.) or type case (upper or lower case), in order to show different population, municipal class, and so on. These differences should be indicated in the legend, normally by displaying an examplar label for each different classification, or possibly by varying the style of the item text itself in the legend.


texts: Python error in QGIS Qfield QFIeldSync plugin configure layers

texts:
I am running QGIS QFieldSync plugin on QGIS 2.14.15. It works wonderfully for packaging my data for QField but I have encountered a python error. When I attempt to package the project for QField using the QFieldSync plugin> Package Project button and click on 'Configure Project Now' if I change any of the layers to 'Offline editing' I get the following python error:

Traceback (most recent call last):
              File "C:/Users/hpethen/.qgis2/python/plugins\qfieldsync\dialogs\project_configuration_dialog.py", line 119, in onAccepted
                self.project.setDirty()
            TypeError: QgsProject.setDirty(bool): not enough arguments

It does not appear to match which layer I attempt to configure to 'Offline Editing'. If I leave them all in 'copy' the programme works correctly.

Have I missed something or is this a bug?


texts: Calculate field values using random integer derived from count of like features

texts:
I have a feature class which contains features classified by type. What I would like to do in my new field is generate a random number sequence for those like values using a count of the like values. Perhaps the example below better illustrates my question:

Field1  Field2
X       random # between 1&4
X       random # between 1&4
X       random # between 1&4
X       random # between 1&4
Y       random # between 1&2
Y       random # between 1&2
Z       random # between 1&3
Z       random # between 1&3
Z       random # between 1&3

I'm thinking I would just use GetCount for the totals and then random.randint(1, GetCount Value) to populate the new field. What I can't seem to figure out is how to have the count reset itself each time it encounters a new value. Thanks for any help you can offer.


texts: QGIS SAGA split lines at points creates redundant lines

texts:
I have 2 projected (ESPG 2039) layers of roads and points (bus stops and junctions). I snapped the points to the roads layer and this is how it looks:

I then try yo use Saga Split lines by points to split the lines but I get redundant lines, e.g. the green and brown lines in the middle of the picture:

Some of the redundant lines are totally independent and some are part of valid lines. I tried using different epsilon values: 0, 0.001, 0.1, 1, 5 and 10.
0.1 gave the least redundant lines, but still many exist.

Is there a way to prevent the creation of these redundant lines?
Sample data (ESPG 2039) can be found here:  https://www.filedropper.com/roadsandpoints


texts: ArcGIS for Windows Phone with Bing maps and shapefiles

texts:
I have a set of shapefiles that I want to display on a map along with satellite imagery.

I currently use ThinkGeo for this sort of thing on Windows.

I'm writing a Windows Phone app that must operate in offline mode and ThinkGeo do not yet have a Windows Phone product.

I think I have 2 choices:

Use ArcGIS for Windows Phone. Use Bing maps imagery when online. Store the shapefiles offline.
Use Bing maps and try to draw the shapefiles onto it or generate tiles from the shape files.

Has anyone gone down this road? Any ideas?


texts:
I have a SpatiaLite database with some multipolygon layers loaded.

What I would like to do is to merge (not dissolve) 2 layers into 1.
Both layers have the same SRID and same field name (they are a result of the same algorithm).

I'm looking for a SpatiaLite funcion that does the same thing of QGIS -> Union.

Is there a simple comand to do that?

My solution, using @Alexandre's suggestions:

CREATE TABLE merged AS
SELECT * FROM layer1
UNION ALL
SELECT * FROM layer2;

SELECT RecoverGeometryColumn('merged', 'Geometry',  32632, 'MULTIPOLYGON', 'XY');


texts:
Is there a simple, free and open source command-line-accessible (other otherwise programmatically-available) viewer that can be called to visualize a GIS dataset?

I know one could do fancy things like view a raster with matplotlib after converting to numpy with Python, or also simply open a full-on GIS with a file loaded by giving it as a command-line argument, but I'm looking for something very lightweight that gives a simple visualization.


texts: Create feature task could not be completed because coordinates or measures are out of bounds?

texts:
I'm trying to create a layer based on a feature class, in ArcMap (I've created a geodatabase from ArcCatalog and linked it to PostgreSQL). The feature class that I've created is registered as versioned. I start editing and draw my polygons and when I finish sketch I get this error :

  The create feature task could not be completed. the coordinates or
  measures are out of bounds.

My coordinate system is : Morocco North (Degrees)
Tolerance : 0.001 meters
XY Resolution : 0.0001 meters

Any thoughts ?


texts: ArcGIS adding opaque circles around coordinates on map

texts:
I have an array of coordinates (lat/long). I want to create a coloured opaque (so I can see inside the circle) 500m radius circle for each one and then display on a map.

Is there an easy way to do this? I am new to ArcGIS.

I am using arcgis in python (version 1.8.0). So far I have a list coordinates that result in the following through map:

df.spatial.plot(map_widget=...)

I want the green coordinates to become opaque circles with a given radius so that I can see the rest of the coordinates inside them.


texts: ArcGIS feature has both a geographic and a projected coordinate system

texts:
As far as I understood a shapefile is either defined as being a Projected Coordinate System or a Geographic Coordinate System. (Seeing as it is just a bunch of numbers and a coordinate system/metadata to define what they mean)

However, in ArcMap, some of my shapefiles have both projected and geographic coordinate systems defined. (See quoted section below).

What's going on here? And how can I define the projected coordinate system for another shapefile that is only currently defined by the geographic coordinate system without completely converting it and losing the geographic coordinate system?

I find it hard to understand why geospatial processing functions don't all have in-built conversion methods. I'm pretty sure when I was using QGIS it just handles all the coordinate system stuff in the background.

The Layer Properties, Source tab:

  Data Type:  Shapefile Feature Class
  Geometry Type:  Polygon.
  Coordinates have Z values:  No .
  Coordinates have measures:  No .

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

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


texts: Shapefile type description showing attribute error

texts:
import arcpy
esc = arcpy.Describe("C:\Users\Drona_Maps\Desktop\lecture16\shapefiles\world_continents.shp")
print esc.ShapeType

Shape file type description showing attribute error like this:
Traceback (most recent call last):
  File "C:/Users/Drona_Maps/Desktop/ARCPY/1st project/CHAPTER 4 exercise part 1.py", line 3, in <module>
    print esc.ShapeType
AttributeError: DescribeData: Method ShapeType does not exist

Process finished with exit code 1


texts: img and trainData do not overlap when using sperclass in R

texts:
I am working on supervised classification using R. I have GeoTiff (icc_m2_clipped.tif) clipped to some extent as below and vector train data (train_ICC.shp). I have changed by first transforming to same projection system, and reprojected to same crs ("+proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0").
I am not winning when I try to classify the image. I get an error "img and trainData do not overlap". I notice that all train polygons during creation lie within the space of the GeoTiff. Where I am I missing it? Please help.
Also note that I am a newbie in R.
library(raster)
library(sp)
library(rgeos)
library(rgdal)
library(RStoolbox)

C:\Users\CIFORRO-96M75\Documents\stack_question
ICC_m2_clip<- brick("ICC_m2_clipped.tif")
train_ICC<- readOGR("ICC_train.shp", "ICC_train", verbose = T)

#Change projection of train to that of tiff

spTransform(train_ICC, CRS(ICC_m2_clip))
proj4string(train_ICC)<- crs(ICC_m2_clip)

 img_classified<- superClass(, trainData = train_ICC, responseCol = "id",nSamples = 1000, trainPartition = 0.8, model = "rf", kfold = 5, mode = "classification", predict = TRUE )

Error in superClass(ICC_m2_clip, trainData = train_ICC, responseCol = "id", img and trainData do not overlap


texts: ArcGIS Pro 2.4.0 Image Segmentation - Higher spectral and spatial detail for high-resolution images?

texts:
I was watching a video regarding image segmentation in ArcGIS Pro 2.4.0. It was an ESRI event (2018), so the presenters were ESRI employees. Here's the video I am referring to: https://www.youtube.com/watch?v=v40abIwGZWA

During the demonstration, they segmented an image using the Segmentation tool under the tab Imagery --> Classification Tools --> Segmentation. The image was a satellite image (I forget exactly which satellite it was from, but I think the spatial resolution was 50 cm/pixel). With this resolution, they said that the default settings for spatial and spectral detail were just fine to use. Below is an image with the default settings:

I am trying to perform image segmentation on very high-resolution drone imagery, with a spatial resolution of 1.1 cm/pixel.

My question is this: Can I have higher spectral and spatial settings for my imagery than the defaults? How would that affect my image classification?


texts: Calculating frequency and area percentage of certain type of landcover in raster data using ArcGIS for Desktop?

texts:
I am a rookie in ArcGIS and I am learning raster.

My problem is that I have a area of different kinds of landcover(actually it's rock type). Now the road is set to NoData, therefore the whole area is divided into several regions/blocks.

Is there a way for me to get know the area percentage of the most frequent rock type within each block?

This is a theoretical problem. My thought is to get the most frequent rock type, I need the help of attribute table-the one with the most count shall be the winner. For percentage I will just need to use the counts of cells of that rock type divided by the number of all the cells of that region.

Is there any way for me to get the total count of the block?
There must be a better way of doing this using certain tool than stupidly count cells and do calculation.


texts: Batch process doesn't work with multiple files

texts:
I'm attempting to determine land surface temperature from Landsat 5 in ArcMap 10.2.2. I created a model to generate land surface temperature. The model works correctly when I run any single landsat scene (date) through the model. When I batch the model with more than one scene I get outputs that are the same to the hundredths place, and which are different than if I ran the scenes individually through the model. The temperature ranges are still reasonable, but they are wrong. If I run ~30 scenes at the same time through the model I get an error about the model trying to access a virtual object. This always causes ArcMap to crash.

I can run the model with each of my scenes one at a time, but I have dozens of scenes I need to process. Is it a quirk of ArcMap that the model runs correctly with one scene and then acts progressively stranger with more scenes added, or is there some setting that I'm missing somewhere?


texts: Single pie chart out of set won't render QGIS

texts: I've joined a CSV with a shapefile to demonstrate boardings/alightings at specific transit stops, and was able to successfully generate a pie chart for all but one of the points via the Diagrams tool. I haven't found another question similar to this one on the forums. Any ideas would be helpful. I've included an image of the map.


texts: How can I move the toolbar if its holder is out of screen?

texts: While I was editing the toolbar, the holder stayed out of the screen. so I can't move it. How can I fix it?


texts: Orthophoto in QGIS to georeference for use in Sketchup 8, or how to match it with Sketchup 8 pictures?

texts:
I have a high resolution orthophoto of an area I want to rebuild in sketchup and therefore want to lign it up (georeference it properly), but all the scaling and everything I tried so far did not help matching it exactly up with the sketchup 8 provided google earth pictures.

I figured out, that it must be due to the fact, that my pictures are in the ETRS89 UTM32N date and the sketchup photos in another, I assume it should be Google Mercator as I researched so far.

(I tried and succeeded matching it using the openlayers plugin in QGIS and the on the fly function, but when I zoom in, everything is out of order again.)
And even if I would not have the zoom in problem, how could I export the layer with the new projection?

I think it´s not possible in QGIS, since you can only save project files.
So is there a nother capable freeware I could use?

Please help me,
any hints very much appreciated.
Thank you for your help,
Os


texts:
I have rasterfiles and I'm using Extract Multi Values to Point but there are some parts of two layers aren't overlap. Normally, these points should get -999 or something like this? But my NoData points getting zero and this could be a major problem because I couldn't understand which one is real zero and which one is NoData.

How can I solve this?


texts: Editing Parcel Fabric in Versioned Geodatabase

texts:
We keep our parcel fabric in a SDE database with an intermediate staging version that multiple users post changes to from their own personal edit versions.

We're able to edit the parcel fabric directly on the intermediate version, but when users try to edit the fabric on their own version, the parcel editing tools remain deactivated. Everyone can successfully edit and post any of the other feature classes in the database.

Is it possible to edit the fabric within these personal versions?

Version tree:

Default

QA

User 1
User 2
User 3


texts: Filter a range number in Query Builder of QGIS 2.8.6

texts:
I need to filter features that have attributes that start within a relative number range. I am using brackets to filter the range of numbers, but it does not work me. This is my query:

CASE
WHEN  "nombre"  like '[0-3]%' THEN 'ROSETAS'
WHEN  "nombre"  like '[4-5]%' THEN 'CTO'
ELSE 'OTRAS'
END


texts: Hard to get my head around, generating average price for houses in each quadrant

texts:
I'm able to determine the average price of houses sold in a particular area (like the squares on image below) dependent upon postcode. Although I can generate the average sale price for a particular postcode, I'm having difficulty utilizing these to the generations to the average house price for each quadrant. As in image.

Assuming I know the average house sale price for each of the letters noted on the image below, how would you generate the average price for each colored quadrant?

Further information:

The information contained within each quadrant (red, orange, green) should only impact that quadrant. (It's these values I hope to utilize to form an average for each separate quadrant)
These quadrants are intended to be eventually presented within Google Maps using circles (which I know how to do)


texts: Setting map extent using offset in ArcGIS API for JavaScript v3.x

texts:
I am trying to set map extent using an offset because there is a panel floating over the left side of the map and sometimes, when I zoom to a feature extent, the feature is partly hidden behind it.

I would like the extent to be always displayed and CENTERED in the map section that is not hidden by the panel

For instance the normal behavior is like below, but some of the feature is hidden behind the panel:

I would like to find a way to center the feature in the map section that is not hidden by the panel like the red rectangle below:

This is what I tried so far without success:

function zoomToExtent(extent) {
  var hiddenMapPartTopLeft = map.toMap(new ScreenPoint(0, 0)),
      hiddenMapPartBottomRight = map.toMap(new ScreenPoint(panel.position.left + panel.position.width, map.container.offsetHeight));

  var dx = new Extent({
    xmin: hiddenMapPartTopLeft.x,
    ymin: hiddenMapPartBottomRight.y,
    xmax: hiddenMapPartBottomRight.x,
    ymax: hiddenMapPartTopLeft.y,
    spatialReference: map.spatialReference
  }).getWidth();

  extent = extent.offset(dx, 0);
  map.setExtent(extent);
}

The map zoom way too far to the right.

Is there a solution for this?


texts: "Print" and "ui.Label" return different strings in console and in the UI in Google Earth Engine

texts:
I am working on a variable label in Google Earth Engine and noticed that the label as printed in the console is not the same as the label shown on the map. Why is that? Code below:

var aoi = ee.Geometry.Polygon([[[-106.75701284923326, 42.31332966408458],
          [-106.75701284923326, 42.003856624775146],
          [-106.31755972423326, 42.06913460511273]]]);

var chirpsSelect = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
            .filterDate('1985-01-01','2018-01-01')
            .filter(ee.Filter.dayOfYear(3,6))
            .filterBounds(aoi)

var chirpsStartImage = ee.String(chirpsSelect.first().get('system:index'))
var dateStartString = ((((chirpsStartImage.slice(6,8)).cat('/')).cat(chirpsStartImage.slice(4,6))).cat('/')).cat(chirpsStartImage.slice(2,4))

var sortChirpsSelect = chirpsSelect.sort('system:index',false)

var chirpsEndImage = ee.String(sortChirpsSelect.first().get('system:index'))
var dateEndString = ((((chirpsEndImage.slice(6,8)).cat('/')).cat(chirpsEndImage.slice(4,6))).cat('/')).cat(chirpsEndImage.slice(2,4))
var dateTitle = ee.String((dateStartString.cat(' to ')).cat(dateEndString))

print(dateTitle)

var dateLabel = ui.Label({
  value: dateTitle
})

var panel = ui.Panel()
panel.add(dateLabel)
Map.add(panel)


texts: How to polygonise a large and detailed raster into a single valid polygon with QGIS?

texts:
I am doing research on the protected areas of the Amazon rainforest and I need to create a single polygon layer that has all the forested non-protected areas in one attribute row. I have tried using different polygonise functions in QGIS with different dissolve functions and fixing geometries, but I always end up with the same final layer which has invalid geometry.

After polygonising the raster and saving the forested areas, I have a multipolygon layer (with thousands of rows in the attribute table) which looks like what I wish to achieve:

However, when I try to dissolve the layer into a single polygon layer (with SAGA Polygon Dissolve (all polygons), but the same happens with GDAL Dissolve), I get these polygons within polygons:

(This is just one area, not the whole layer, showing the areas that appear after trying to dissolve the multipolygon layer)

I've tried to use "Fix geometries", but that only removed some of the unwanted areas (the fixed image is included here).

I have tried to fix this problem for more than a week now and am running out of options to try.

Edit:
I am running 3.2.0-Bonn on Windows 10. My GDAL/OGR version is 2.2.4 and I use WGS84 as my CRS.


texts: Using SQL in QGIS field calculator to populate new column?

texts:
Courtesy of Geofabrik I have a "places" point shapefile which includes in the attribute table a column named "type". As can be seen in the attached image the column contains fields such as "city", "town" and "village".

For labelling purposes, how do I populate the "font" column with integers such that when "type" = 'village' then "font" will be set to (say) 7, or when "type" = 'city' then "font" will be set to (say) 12, and so on?


texts: Fixing Script Error Invalid Pointer from ArcGIS Geoprocessing Tools?

texts:
I am just learning ArcGIS 10.2 and it was working well until I tried to use the Export to a Geodatabase (single) option for a dbs table I made. Since then, it seems I now get the same Script Error for import/export functions. However, I just tried to export to Geodatabase (multiple) and error did not pop up.
I spoke with tech support and they could not help me so far. I really do not know much about GIS so it is hard for me to fix this problem and cannot proceed with some things.


texts:
How can I do spatial selections or filters in R? I can do this in PostGIS using (e.g.):

SELECT * FROM mypoint, mypoly
WHERE ST_Within(mypoint.geom, mypoly.geom);

and it is similarly a piece of cake in most other GIS software. So I'd expect something in R to look like:

mypoint[within(mypoint, mypoly),]

but I can't find any leads anywhere. I'm looking for something that uses the common spatial predicate terminology used in PostGIS/GEOS/JTS/Shapely (within, contains, intersects, overlaps, touches, etc.; DWithin would be bonus). Does this functionality or library exist for R? If not, is there something that is similar?

Note: In R, I use the sp class for my objects (SpatialPointsDataFrame, etc).


texts: how to install plugin through the python console in QGIS

texts:
I have been trying to load/install plugins via the python console in QGIS. After going through the documentation, I'm sure this is the way. I want to check if a plug-in is already installed, if not install it.

The lines:

from qgis import utils

try:
    utils.isPluginLoaded('zoomtocoordinates')
except ValueError:
    qgis.core.QgisInterface.addPluginToMenu('zoomtocoordinates')

The error returned:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'module' object has no attribute 'QgisInterface'


texts: Restoring the exact original situation of a spatial bookmark?

texts: After saving a view with a spatial bookmark, if one later wants to get back to this exact same view, it appears that you have to have the correct height/width ratio in the view panel. Is there a way to save and reestablish that?


texts: How to call QgisInterface class in a plugin

texts:
I'm trying to access the QgsLayerTreeView class through the QgisInterface class in the python plugin I have created. Ultimately I want to be able to return the selected node (group or layer) from the QgsLayerTree.
The issue I'm running into is I don't know how to call it. Most information suggests to use iface as a way to access the interface but it's not recognised when called from the plugin (works fine from the python console in QGIS).
Is there any other way to create an instance of the QgisInterface class? Or is there a way I can import the iface class to work in my plugin?


texts: Adding features to feature class that is initially empty using ArcPy?

texts:
I have built a loop that creates each time a series of features (lines). And every time that new features are created I want to save them into a feature class (always the same one).

I believe I need to use Append_management to do this. The problem is that the target feature class will initially be empty and therefore the first time the script is trying to append things to it I get an error stating that it doesn't exist.

#Import geoprocessing.
import arcpy
import os

#sys.path[0] gives us the folder where the script is. Then we have to go up one level in the folder structure
#(first we build the path to the parent directory and then we need to normalise it
working_folder = os.path.normpath(os.path.join(os.path.dirname(__file__), '..'))

#Input variables.
fc_Opt_All = os.path.join(working_folder, 'DATA\Optimity_Equipment.gdb', 'Optimity_All')

#Building all existing Wibre links
#Selecting features from Optimity_All where 'LINK' field is not blank and saving them into a feature class
arcpy.MakeFeatureLayer_management(fc_Opt_All, 'fl_Opt_All')
query = "LINK NOT LIKE ''"
arcpy.SelectLayerByAttribute_management('fl_Opt_All', 'NEW_SELECTION', query)
fc_airB = os.path.join(working_folder, 'DATA\Optimity_Temp.gdb', 'airB')
arcpy.CopyFeatures_management ('fl_Opt_All', fc_airB)

#Building all links between airA (indicated in field 'LINK') and airB
fc_links = os.path.join(working_folder, 'DATA\Optimity_Temp.gdb', 'links')
field_names = ['LOCATION', 'ESTATE', 'ADDRESS', 'OFFSETA', 'OFFSETB', 'Z', 'OFFSETC', 'AZIMUTH1', 'AZIMUTH2', 'VERT1', 'VERT2', 'RADIUS1', 'RADIUS2', 'STATUS', 'LINK']
with arcpy.da.UpdateCursor(fc_airB, field_names) as cursor_1:
    for row_1 in cursor_1:

        arcpy.MakeFeatureLayer_management(fc_airB, 'fl_airB')

        #Selecting airB
        query = "LOCATION = '{}'".format(row_1[0])
        arcpy.SelectLayerByAttribute_management('fl_airB', 'NEW_SELECTION', query)
        arcpy.AddMessage(query)

        #Selecting airA
        arcpy.MakeFeatureLayer_management(fc_Opt_All, 'fl_Opt_All')
        query = "LOCATION like '{}%'".format(row_1[14])
        arcpy.SelectLayerByAttribute_management('fl_Opt_All', 'NEW_SELECTION', query)
        arcpy.AddMessage(query)

        #Building the line between the pair
        fc_link = os.path.join(working_folder, 'DATA\Optimity_Temp.gdb', 'link')
        arcpy.ddd.ConstructSightLines ('fl_airB', 'fl_Opt_All', fc_link, 'OFFSETC', 'OFFSETC')

        arcpy.Append_management(fc_link, fc_links)

The error message I get is:

Traceback (most recent call last):
  File "C:\ESRI\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\Affected_LOS.py", line 43, in <module>
    arcpy.Append_management(fc_link, fc_links)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 3984, in Append
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 3981, in Append
    retval = convertArcObjectToPythonObject(gp.Append_management(*gp_fixargs((inputs, target, schema_type, field_mapping, subtype), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
 ERROR 000732: Target Dataset: Dataset C:\ESRI\ArcGIS_Pro_Projects\Optimity_Alf\DATA\Optimity_Temp.gdb\links does not exist or is not supported
Failed to execute (Append).
 Failed to execute (AffectedLOS).

Any idea on how to get this working?


texts: ST_intersection returns erroneous polygons

texts:
I am trying to perform a clip by attribute in PostGIS between two shapefile containing 3000 polygons.
Both shapefiles are free of errors (checked with topology checker), both shapefile are using EPSG:3857.

The following query works well for 99.9% of polygons, but for some it creates incomplete polygons or misses some when loaded in QGIS. The polygons with errors do have the same name used as a clipping attribute. I really don't understand why some polygons fail to be generated compared to others.

create table poly_20m_checked (gid serial PRIMARY KEY, polyID text, geometry geometry);
    insert into poly_20m_checked(polyID, geometry)
    select
    a.name as polyID,
    ST_Intersection(a.geom,b.geom) as geometry
    FROM donut20m a, symm20m b
    WHERE a.name=b.name
    Group by polyID, geometry

donut20m:

symm20m:

result (one polygon is missing):

Sample files that contains 4 polygons, among which one that returns an error (moved to make location anonymous):
https://drive.google.com/open?id=0B7mqg7Gq1NGVVVVsbTZqMlVhZnc

The reason I am doing it in PostGIS is that I don't think there's a way to clip polygons in 2 shapefiles by using attribute in QGIS. What could be causing errors?
I am very open to suggestion to get around the problem. I thought about spliting polygons by ID and then clip each pair of shapefiles but haven't really managed to code it.

EDIT: download for sample polygons.
EDIT2: Tried to simplify geometries by 1m, which makes thing a bit better, but still bugy


texts: Implementing MultiValue Parameter with Multiple Columns with Python and ArcGIS Desktop?

texts:
I am working on an ArcGIS (9.3) Python script. One of the inputs is a list of features. For each feature, I need to select a corresponding field. Setting the parameter to MultiValue gives a nice table in the UI for all the features being processed. How do I add another column to the MultiValue table to select a field for each feature?

Cross-posted: http://forums.arcgis.com/threads/30907-MultiValue-Parameter-with-Multiple-Columns

Esri's Answer from the ArcGIS forum:

  Unfortunately, the solution for this was implemented just recently for
  10.1 with ValueTable inputs for script tools. You will need to have two multivalue parameters with some custom validation to make this
  work for 9.3.1.


texts: Can't join coordinates to polygons represented by black icons in RStudio viewer

texts:
I'm using {sf} to import a GeoJSON file. What does the black icon represent? I'm trying to join coordinates to this object via st_join(ken_con_sf, join = st_within) and it seems to fail where coordinates fall within the polygons represented by these icons.

library(sf)
ken_con_sf <-
    read_sf("https://raw.githubusercontent.com/mikelmaron/kenya-election-data/master/data/constituencies.geojson") %>%
    st_transform(4326)

R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] maps_3.3.0           sp_1.4-2             sf_0.9-6
 [4] ggrepel_0.8.2        zoo_1.8-8            leaflet_2.0.3
 [7] showtext_0.9         showtextdb_3.0       sysfonts_0.8.1
[10] lubridate_1.7.9      shinyWidgets_0.5.3   rmarkdown_2.5
[13] shiny_1.5.0          flexdashboard_0.5.2  dplyr.teradata_0.3.2
[16] niviR_0.1.0          dbplyr_1.4.4         forcats_0.5.0
[19] stringr_1.4.0        dplyr_1.0.2          purrr_0.3.4
[22] readr_1.4.0          tidyr_1.1.2          tibble_3.0.4
[25] ggplot2_3.3.2        tidyverse_1.3.0      RSQLite_2.2.1
[28] RPostgres_1.2.0

loaded via a namespace (and not attached):
 [1] fs_1.5.0           bit64_4.0.5        httr_1.4.2
 [4] tools_4.0.2        backports_1.1.10   utf8_1.1.4
 [7] R6_2.5.0           KernSmooth_2.23-17 DBI_1.1.0
[10] colorspace_1.4-1   withr_2.3.0        tidyselect_1.1.0
[13] bit_4.0.4          compiler_4.0.2     cli_2.1.0
[16] rvest_0.3.6        xml2_1.3.2         bookdown_0.21
[19] scales_1.1.1       classInt_0.4-3     odbc_1.2.3
[22] digest_0.6.27      pkgconfig_2.0.3    htmltools_0.5.0
[25] fastmap_1.0.1      htmlwidgets_1.5.1  rlang_0.4.8
[28] readxl_1.3.1       keyring_1.1.0      rstudioapi_0.11
[31] generics_0.1.0     jsonlite_1.7.1     crosstalk_1.1.0.1
[34] magrittr_1.5       Rcpp_1.0.5         munsell_0.5.0
[37] fansi_0.4.1        lifecycle_0.2.0    stringi_1.5.3
[40] yaml_2.2.1         grid_4.0.2         blob_1.2.1
[43] promises_1.1.1     crayon_1.3.4       lattice_0.20-41
[46] haven_2.3.1        hms_0.5.3          knitr_1.30
[49] pillar_1.4.6       reprex_0.3.0       glue_1.4.2
[52] evaluate_0.14      modelr_0.1.8       vctrs_0.3.4
[55] httpuv_1.5.4       cellranger_1.1.0   gtable_0.3.0
[58] assertthat_0.2.1   xfun_0.19          mime_0.9
[61] xtable_1.8-4       broom_0.7.0        e1071_1.7-3
[64] later_1.1.0.1      rsconnect_0.8.16   class_7.3-17
[67] memoise_1.1.0      units_0.6-7        ellipsis_0.3.1


texts:
I am writing a plugin which shows a graph in a qGraphicsView based on some selected points. I want to be able to read out points clicking in the graph, so I've written a mosePressEvent eventhandler but it does never fire. To test it, I just did

def mousePressEvent(self, event):
      self.iface.messageBar().pushMessage(
                "Success", "Pressed",
                level=Qgis.Success, duration=3)

But no matter where I click on the widget, no message is shown. I do neither get any error messages. Do I need to import something or do something else special to make it work?

(The plugin is based on a template by the Plugin builder plugin, maybe I need to remove something)

The entire file is in https://github.com/sickel/qgisSpectre/blob/1beff9842ededc425d3200d566d11dde9c1d3bb7/qgisSpectre.py


texts: Looking for a FREE Cloud based 'Web Mapping Stack' Host

texts:
Any pointers to free cloud-based geoweb hosts? Use case: To display a basic interactive web map in Leaflet/ OpenLayers with data hosted in PostGIS served by GeoServer. I am looking at something very close to the OpenGeoSuite since I am playing with that at the moment already. [This is basically a teach-yourself webmapping exercise and there are zero personally hosted server resources on my sleeves].

I've found capnix to host the website but now need the 'Geo'-portion. CartoDB's 5MB data limit is too small for my project needs and OpenGeo's Cloud Edition is too short (3 days!). Any ideas?


texts:
Is there a way to run st_split from within QGIS?  I am looking to split a polygon with a line feature.  My data is stored in a PostGIS database, and I would like to perform this operation from within QGIS.  I am aware of PG_Query for QGIS, but st_split is not a tool within the plugin.

Thanks...


texts: Loading point shapefiles from Trimble GPS in QGIS

texts:
I had some field work with Trimble gps. I located some positions and I exported them as a shapefile in Trimbles pathfinder software, but when I load vector data point shapefile that represent positions, they are off the map.

My map is in WGS 84 CRS; other shapefiles are ok, but these positions are wrong. Any help?
Here are the screenshots