texts: Add a Google Earth Enterprise globe (from geserver) to QGIS

texts:
I am managing Google Earth Enterprise, where I can add data to Fusion, push it to geserver and then publish it to Google Earth Enterprise clients on people's computers. But the people need to add a database by typing in a URL containing e.g. geserver.it.firm.com/globe1/ and no other credentials are needed. I'm assuming credentials are indirectly processed because they need to login to their computers.

I'm wondering if it's possible to add this data into QGIS?
I've tried adding Layers using WMS --> New Connection -->.. for name i just put 'Globe1' and for URL I put http://geserver.it.firm.com/globe1 but it says "Failed to download capabilities"

Is it possible to add this data into QGIS somehow?


texts: How can I use relative paths in print composer for images?

texts:
I'm sharing an entire map project with a client so they can eventually print off an atlas. I've managed to get all the layers and icons to work correctly with relative paths, but the atlas contains images. These don't appear when the client opens the project. Meanwhile, if I truncate the image path to a relative one, the images won't load in my local project.

This seems to be an updated version of this question, and I would also like to avoid modifying startup.py. I notice that I could create an expression to use for the image path, is there some combination of functions that would work for this?


texts: How to rank according to a distance in a shapefile using R?

texts:
I would like some help on the theory behind this. I have done this in GRASS easily but doing it in R is not the same.

Here's the table with the info, there are more tables similar to that but refer to different types. If I get this then I will be able to do the others since it the logic is alike.

%%UPDATE%%

I want the region (buffer of 500) to be ranked with 0. the 500-1000 (buffer of 1000 from the wetlands), 5 and the 1000(rest of land-more than 1000m) with 10

   type      distances    ranking
               <500          0
 wetlands      500-1000      5
               >1000         10

I have imported the shapefiles in R but I don't know what to do with the table?

EDIT(this is what i did so far)

coastline<-readShapeSpatial("shp")
wetlands<-readShapeSpatial("shp")
wetbuf<-gBuffer(wetlands,width=500,byid=T)
wetbuff<-gBuffer(wetlands,width=1000,byid=T)
wetbuf$val=0
wetbuff$val=5
coastline$val=10


texts: Numbering polygons according to their spatial relationships?

texts:
I have, for example, this layer in postgresql (polygons):

table : poly (objectid: integer, num: integer, shape: geometry)

I want to enumerate the polygons from left to right and top to bottom like this:

I am working with Vb.net, ArcEngine, PostgreSQL, and ArcSDE.

Solution that I proposed

I think to create grid lines (red layer; see picture) and iterate it, then verify intersection and finally enumerate the polygons. Will this work?


texts: Just completed install of QGIS (2.0.1 Dufour) on a Macbook Pro running on Mavericks via Kyngchaos download page. However, when downloading Matplotlib for the install, while everything ended up the correct python site package (2.7), the matplotlib wound up in a folder within the 2.7 site package called matplotlib-override. Does anyone have an explanation? When in the 'override folder QGIS could not find Matplotlib, but once removed and bumped up a level to the site package (2.7) folder QGIS installed cleanly (still running it through its paces for functionality).


texts: What type of object is required for IViewRefresh.RefreshItem "Data" parameter?

texts:
Esri's help for this method specifies this as LPUNKNOWN, which is rather unhelpful since I don't speak COM.  I've tried an envelope, a geometry, and a graphicsContainer; each time receiving "System.ArgumentException: Value does not fall within the expected range.".  Can't find any meaningful examples.  Here's my MCVE . . .

IMap newDataFrame = new Map
{
    Name = "New Data Frame",
    SpatialReference = CurrentMap.SpatialReference
};
var viewRefresh = (IViewRefresh) newDataFrame;
var graphicsContainer = (IGraphicsContainer) newDataFrame;
var pointMaker = IdioticThingWithWhichYouNeedNotConcernYourself();
CurrentMxDocument.Maps.Add(newDataFrame);
CurrentMxDocument.UpdateContents();
...
for (int i = 0; i < 5000; i++)
{
    IMarkerElement element = pointMaker.MakePoint(out var envelope);
    graphicsContainer.AddElement((IElement)element, 0);
    viewRefresh.RefreshItem(graphicsContainer);
    // Also tried viewRefresh.RefreshItem(envelope);
}

I figure it's questionable (at best) to refresh even a tiny portion of a map 5000 times, but this is just tinkering so please humor me on that point. Regardless, the exception noted above is thrown the first time through the loop.


texts: Error "no method for coercing this S4 class to a vector"

texts:
I am trying to reproduce the example posted on

Extracting intersection areas in R

Although, when I try to create the SpatialPolygon "p2", R returns the error:
"Error in as.vector(x) : no method for coercing this S4 class to a vector".

Any suggestions why is it happening? How can I fix it?


texts: Having problem with pgsql2shp function - get postgis_version() function not found error

texts:
I try to use pgsql2shp postgis function to export a table from another schema to .shp format.
Base on this link, I run following code
pgsql2shp -f "c:/path/to/export/file" -h localhost -u <username> -P <password> atlas country.parcel

I created a country schema and then added postgis extension to it but I get this error
ERROR:  function postgis_version() does not exist
LINE 1: SELECT postgis_version()
           ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I'm pretty sure that I've added postgis extension to it. I followed some solutions like this but they didn't work. I think it is not related to search_path as it's mentioned in above solution.
What is the problem?


texts: ggmap: create circle symbol where radius represents distance (miles or km)

texts:
I am trying to add a circle around a given latitude and longitude point where the circle size (radius) reflects distance in miles or km.

I managed to come up with this example using ggmap, where I can add a circle. but what do I need to do to have it represent distance. I would ultimately like the bounds of the circle to represent 20 miles from the center.

Code below:

library(ggmap)
library(ggplot2)

d <- data.frame(lat=c(33.79245),
                lon=c(-84.32130))

emory <- get_map("Atlanta,Georgia", zoom=12)

p <- ggmap(emory)
p <- p + geom_point(data=d, aes(x=lon, y=lat),
                    color="red", size=20, alpha=0.5)
p


texts: Evaluating training samples for supervised image classification in ArcGIS Pro?

texts:
Its my first time using ArcGIS Pro and I have started doing a supervised classification.

In most software you have some tools such as histograms, scatterplots and/or statistics to evaluate training samples but I couldn't find any of these tools in the ArcGIS Pro image classification options.

Are they hidden in another menu or did they really didn't include these tools?

I could not find any information regarding this in the help documentation, tutorials or forums.


texts: Selecting which polygon's attributes to keep after dissolve

texts:
I have two polygons, after doing a dissolve I can't choose which polygon's attributes the new dissolved polygon will have. I tried to find a pattern by checking if the bigger's polygon attributes are kept but it does not appear to be the case. The selection seems random. What do you think?

The QGIS version is 3.4.11-Madeira

I suppose the reason it keeps the 63 is that it is higher in the attribute table as a row? Therefore it is normal that it keeps the 63?

EDIT

Here in the one layer the dissolve keeps the number that is in the higher position as a row in the attribute table which implies that prevalence is the indicator of which id it will eventually have. But the second layer that has been sorted keeps the same id as the other layer after the dissolve which means that the position in the attribute table may not be of importance in this case.


texts:
I am receiving an error ("Error executing identify") when I try to perform an identify operation on map running on ArcGIS Server 9.3.  All of the inputs appear to be correct when I check them in the ArcGIS Server Directory.

The map I am using has over 90 layers.  Could that possibly be causing an issue with the identify operation?


texts: How can I apply a reverse QgsCoordinateTransform on a QgsGeometry in pyqgis?

texts:
I have defined xform as a QgsCoordinateTransform to project geographic coordinates.
I can use unprojPoint = xform.transform(point, QgsCoordinateTransform.ReverseTransform) to do the reverse operation, but how can I apply this to a QgsGeometry ?

I can do line.transform(xform) allright but this won't take a direction argument.
I would expect something along the lines of line.transform(xform, QgsCoordinateTransform.ReverseTransform) would work but no luck...

Do I have to setup a new QgsCoordinateTransform to do this ?

Note: I am using Qgis 3.0


texts:
I have a shp containing polygons with ELEVATION attribute field. I want to import it to Autocad and display it as 3D polygons based on the ELEVATION field without using SOLID or SURFACE in acad .
I've seen that "Spatail Manager" (http://www.spatialmanager.com/spm-forautocad/) gives a solution but with cost of payment.
is it possible to do it with Autocad only ?

This video shows what I'm looking for: https://www.youtube.com/watch?v=dzSNkOBtQnU&t=2s

-I'm using AutocadMap3D 2018.


texts:
I am looking for a tool from the GDAL that is able to reclassify a raster containing multiple discrete values. So far, I have found Reclassify rasters using GDAL and Python; however it seems that only single values can be handled here.

Is there something like a native tool from the suite?


texts: What does GDAL's function addo stand for?

texts: This is sort of a superficial question. I can't seem to find the meaning behind the naming convention of gdaladdo. My guess is that addo stands for "add overview"? Does anyone else know or have a better guess?


texts: Crop a raster image with the coordinates in a CSV file in Python

texts: I am trying to crop a satellite image with geographical coordinate from  a CVS file. I have never cropped an image using coordinates from CSV file. I tried searching a lot but I can only find cropping using a shapefile. Can someone guide me or give me a link of tutorial  or something where they have cropped using a CSV file. I really don't know where to look.


texts:
I downloaded the map of my city already in shapefile into postgis http://metro.teczno.com/ to load in order to consult routes using pgRouting.

Iss happening the error: http://metro.teczno.com/

Importing with configuration: port-alegre.osm online, public, geom, C: \ Cassia \ Shape POA-OSM \ port-alegre.osm-line, mode = c, dump = 1, simple = 0, geography = 0, index = 1, shape = 1, srid = 4326
C: \ Cassia \ Shape POA-OSM \ port-alegre.osm online: DBF file (. Dbf) can not be opened.
Shapefile import failed.

I do not know what the geometry column in esquecificar postgisloader, can someone help me?


texts: configurable map viewer multiple map services

texts:
I have multiple map layers (right now about a dozen but that will likely stretch to several dozen) in my CMV. I have introduced several WMS services to operationalLayers that I have published to a local server. Only the "Parcels" service will show up. I have tried different combinations of several layers but I can never get more than one layer to display at a time. Does anybody know why multiple services will not display in my CMV?

This is what the section in my viewer.js looks like:

operationalLayers: [{
        type: 'dynamic',
        url: 'http://arcgis.mydomain.com/arcgis/rest/services/Floodplain/Floodplain/MapServer',
        title: '100yr Floodplain',
        options: {
            id: 'Floodplain',
            opacity: 1.0,
            visible: true,
            outFields: ['*']
        },
        identifyLayerInfos: {
            layerIds: [0,1,2,3,4]
        },
        legendLayerInfos: {
            layerInfo: {
                hideLayers: [21]
            }
        },
        type: 'dynamic',
        url: 'http://arcgis.mydomain.com/arcgis/rest/services/Parcels/Parcels/MapServer',
        title: 'NY Parcels',
        options: {
            id: 'Parcels',
            opacity: 1.0,
            visible: true,
            outFields: ['*']
        },
        identifyLayerInfos: {
            layerIds: [0,1,2,3,4]
        },
        legendLayerInfos: {
            layerInfo: {
                hideLayers: [21]
            }
        }
}],


texts: How to divide attribute table in map composer between 2 or more pages?

texts: Does anyone know how to divide attribute table between 2 or more pages? In exported PDF it's unfortunately clipped.


texts:
I have a query filter which works if pass a new instance of QueryFilterClass without setting any properties.  But it breaks if I set the WhereClause.  I get an error
0x80040358

queryFilter.WhereClause = "HouseNum = '4501'"; //this is my WhereClause which will make table.Search(queryFilter, true) break

This is my code which works:

IQueryFilter queryFilter = new QueryFilterClass();
//queryFilter.WhereClause = "HouseNum = '4501'";

IQueryFilterDefinition queryFilterDef = (IQueryFilterDefinition)queryFilter;
FileInfo file = new FileInfo(dbfPath);
String dbfDirectoryPath = file.Directory.FullName;
IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactory();
Workspace workspace = workspaceFactory.OpenFromFile(dbfDirectoryPath, 0) as Workspace;
IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;
String dbfTable = file.Name;
ITable table = featureWorkspace.OpenTable(dbfTable);

int shapeIndex = table.FindField("Shape");

try
{
using (ComReleaser comReleaser = new ComReleaser())
{
ICursor cursor = table.Search(queryFilter, true);
comReleaser.ManageLifetime(cursor);
IRow row = null;
while ((row = cursor.NextRow()) != null)
{
object o = row.get_Value(shapeIndex);
IPoint p = o as IPoint;
p.Project(gcs);
double doubleX = p.X;
double doubleY = p.Y;
}
}
}

How can I fix my WhereClause?  Also does QueryFilter support the SQL LIKE Condition?


texts:
Does anyone have ideas or know how one might be able to use CartoDB to create a map that one can toggle the individual features within a KML file, as currently happens on this map?

I want my users to be able to toggle individual routes on and off, and Google is removing the feature that allows the KML/KMZ files to be overlaid on Google Maps in February.

I already have a barebones setup based off some other research I've done, but was hoping that maybe there is a more robust solution out there that is closer to what I have right now.


texts:
I work with this polygon and try to find a way to create a new polygon from the empty gap (not distinguish it) in the middle of the polygon.

I succeed to do it with a "Difference" tool but used with another layer on the polygon and then using the difference tool.

I can also digitize the shape of the hole and create a new layer I seeking an automated way to do it and not doing in manually.

I read:

"Clip" hole into Shapefile
Create Polygon from holes in polygon ArcGIS
How to fill in a ring polygon in QGIS 2.2?

Is it possible to do it with a direct tool in QGIS 2.18.1, without the aid of another layer?


texts: Use relation to apply same colour to children feature

texts: I'm looking for a function with field calculator in a symbology context to apply the colour of parents features to the children with transparency.


texts: Create polygon from lines in PostGIS : Problem

texts:
I create a new post to clarify more. My former answer was not adapted here : (Create polygon from lines in PostGIS
I apologize for my english.

I have a problem to build polygons.

And I run a query :

CREATE TABLE "CTRL_SONDE_XXXX"."067_georezo"
AS
  WITH t AS (
    SELECT
      p.id_triangle,
      array_agg(l.geom)||st_exteriorring(p.geom) AS t
    FROM "CTRL_SONDE_XXXX"."067_TRIANGLE_non_intersect_Delau" p,
      "CTRL_SONDE_XXXX"."067_line_intertriangle_Delau_step1_sel" l
    WHERE ST_Intersects(
      ST_SnapToGrid(p.geom,0.0000001),
      ST_SnapToGrid(l.geom,0.0000001)
    )
    GROUP BY p.id_triangle,p.geom
  ),
  u AS (
    SELECT id_triangle, st_dump(st_node(st_collect(t))) d
    FROM t
  )
  SELECT (st_dump(st_polygonize((d).geom ORDER BY (d).path))).geom
  FROM u
  GROUP BY id_triangle;

but the result is wrong :

Do you have any idea?


texts: Join multiple Points to Polygon and concatenate joined point values of column into Polygon column

texts:
My data:

Points:
    ID HN Street
    1  5  Examplestreet
    2  6  Examplestreet
    3  2  Otherstreet
    4  2  Anotherstreet
Polygons:
    Name
    Firstpolygon
    Secondpolygon
    Otherpolygon

So I want to performe a one(polygon) to many (points) spatial join and concatenate the joined valus from the points HN column to a new column "all_HN" in the Polygonsfile.

So the result should look like this:

Name          all_HN
Firstpolygon  5,6   -> if the points ID1 and ID2 lie within the same polygon (`"Firstpolygon"`)
Secondpolygon 2     -> point ID3 within the `"Secondpolygon"`
Otherpolygon  NULL  -> no point within "Otherpolygon"
and so on

I thought of solving this with geopandas. (I have about 200.000 Points and 100.000 polygons)
Using the following code:

from geopandas import gpd

points = gpd.GeoDataFrame.from_file('MyPointsFile) # or geojson etc
polys = gpd.GeoDataFrame.from_file('MyPolygonsFile.shp')
pointInPoly = gpd.sjoin(points, polys, op='within',how='inner')

Now I thought using something like:

pointInPoly.groupby('index_right')['HN_left'].sum()

But instead of sum() the right command to concantinate all values of the matching points into a new column "all_HN".

Anyone, any idea how to solve this problem?Its also ok using another package and not geopandas.
Since geopandas is based on pandas an pandas solution should also work.
It would also be enough to just append all columns of the matches to the corresponding Polygons Tables.


texts:
I am using qgis 2.2 (I have also tried this in qgis 2.0 and arcgis 10 with the same result) (I have some basic understanding of netcdf, but I am probably missing a lot of details...)

I have manged to read in some netcdf-files and have them displayed correctly. Some other fails miserably, as the coordinate infomation is ignored and the file is displayed south- and eastwards from 0,0 ... Is this a known problem? I have been looking at the files with ncdump, but have not been able to spot any differences that should make one of them work and the other to fail. When looking at the metadata in qgis or arcgis, it is clear that a lot is missing for those files that do not display correctly.

Edit:
When converting the files from netcdf4 to netcdf3 they worked fine. is this a known issue in arcgis / qgis? (using arcgis 10, it seems to use gdal, but I cannot find which version, qgis 2.2 from osgeo4w using gdal 1.10)

Edit2:
The original files were so large it was a bit difficult to share them. Now I have made an example file that behaves the same way: http://sickel.net/div/test.nc - it is created from http://sickel.net/div/test.cdl - it should, as far as I do understand, cover an area from -160 to -15 E, 20 - 60 N ie most of US main land and neary to UK in the atlantic ocean, but it ends up from 0,0 and south and eastwards in qgis, in arcmap 10.0, it ends up as a point on -160,20 ...


texts: is it possible to combine multiple URLs for one ol.layer.tile in openLayers 3?

texts:
I have two Image services that I would like to combine into one layer; one has higher resolution imagery that draws only at a larger scale, and other has lower resolution imagery that draws at all scales. I'd like to combine them so that when you are zoomed in, the detailed imagery draws on top of the lower res imagery. I tried to load in the 2 urls like this:

var USGSimagery = new ol.layer.Tile({
    myattribute:'USGSimagery',
    source: new ol.source.TileWMS(({
    urls: ['http://raster.nationalmap.gov/arcgis/services/Orthoimagery/USGS_EROS_Ortho_SCALE/ImageServer/WMSServer','http://basemap.nationalmap.gov/arcgis/services/USGSImageryOnly/MapServer/WmsServer'],
    params: {
                  'LAYERS': 0
                        }
       }))
});

But when I load the layer in the map - it draws funny, like in alternating strips where it shows one service and then the other. So this obviously isn't the way to do it. Does anyone know if it's possible to combine two services into one layer with one drawing on top of the other?


texts: Exporting GeoJSON from QGIS to Mapbox Studio

texts:
I am willing to create GeoJSON layer with QGIS, download it and upload it to Mapbox Studio. When I upload the GeoJSON file to Mapbox Studio, I get the following error:

Input failed. Datasets don't support features outside longitude
+/-180, latitude +/-90.

This is how I create the GeoJSON File:

Create new Point shape layer
Add Features
Set CRS to Pseudo Mercator EPSG:3857
Save Layer as GeoJSON

What am I doing wrong here?


texts: Doing spatial join in QGIS between land parcels and zoning map, how can land use with biggest overlap give its attribute to parcel?

texts: The parcel might have some small overlaps (or intersections?) from other land uses if the alignment of vertices is not perfect, I want to make sure the land use covering the largest part of the parcel is the one taken into account by QGIS. Both layers will be polygons. Not sure if I should choose intersect or overlap in this case also.


texts: I'm new at this, a novice at best. I'm using QGIS 1.7.3. I simply want to select a large number of points on a single polygon and delete them. Is this possible?


texts: Is there a name for map projections that produce two disks?

texts:
I've seen lots of world maps in the shape of two disks that touch, each one a separate projection of one hemisphere onto a disk -- usually the Western Hemisphere and the Eastern Hemisphere -- so together the full map covers the entire globe.

Is there a standard name for this kind of map projection onto 2 circular disks?

from Wikipedia: Rumold Mercator

from Wikipedia commons

Ideally I'd like a general name for all "2-disk" arrangements like this,
but I'd be happy with separate specific names for
"2 disks that use a Wulff net" and
"2 disks that use a Schmidt net" and
"2 disks that use a Mollweide map".


texts:
I've written a short python script for use in QGIS (currently using 2.8.1) to set the CRS for all loaded layers to OSGB 1936 / BNG. The script correctly sets the CRS but I have to manually scroll in/out to refresh the screen after running it - not a big issue though slightly annoying and possibly easy to fix!

I've tried using the code snippet below to refresh but I still have to manually scroll in/out.

qgis.utils.iface.mapCanvas().refresh()

Is there a line of code to force a full refresh/re-render of all loaded data?


texts: How do I compute the forest area for different regions?

texts:
I am currently struggling to compute the forest area of some regions. I am an absolute beginner with Google Earth Engine.
I tried it 2 different ways, mostly by copying and editing code I found.

Count the number of pixels with a tree cover greater than 10% (using the Hansen data)
I thought that, I then could simply multiply by the pixel size and get the area.
However, I am not sure if I can do it like that. Nonetheless, here is my code:

    var gfc2017 = ee.Image('UMD/hansen/global_forest_change_2017_v1_5');
        var treeCover = gfc2017.select(['treecover2000']).unitScale(10, 100);

        Map.addLayer(
            treeCover,
           {palette: ['000000', '00FF00'] });

    // import the District layers as a fusion table
    var districts = ee.FeatureCollection('ft:13cHA4yautjvngKSaPywzP9Li9FVtlkJ7Dg64DQMG');
    //add districts to the map and center map
    Map.addLayer(districts);
    Map.centerObject(districts);

    // Count the pixels with tree cover over 10 percent
    var stats = treeCover.reduceRegions({
    collection: districts,
     reducer: ee.Reducer.count(),
     scale: 30
    });
    //print('count of pixels representing forest in each circle:', stats)

    print('pixels representing forest: ', stats);

The numbers I am getting are however, way to high. Can anybody please help me?

My second approach was to try using a map that has pixels of forest area, non forest area and water of the year 2010 (ID:JAXA/ALOS/PALSAR/YEARLY/FNF/2010)

I tried to sum the pixels, however, I am not able to just look at forest pixels. I guess it is pretty easy to do. But reading here and googling a lot has not helped yet.

My code (without the selection of only forest pixels):

 // Load the image from the archive
    var farea = ee.Image('JAXA/ALOS/PALSAR/YEARLY/FNF/2010')

    var treeCover = farea

    Map.addLayer(treeCover)

// import the District layers as a fusion table
var districts = ee.FeatureCollection('ft:13cHA4yautjvngKSaPywzP9Li9FVtlkJ7Dg64DQMG');
//add districts to the map and center map
Map.addLayer(districts);
Map.centerObject(districts);

// Sum the values of forest pixels in districts.
var areaImage = farea.multiply(ee.Image.pixelArea());

// Count forest pixels in districts
var stats = areaImage.reduceRegions({
collection: districts,
 reducer: ee.Reducer.count(),
 scale: 30
});

Would be really great if someone could help me out.


texts:
I am fairly new to Arcgis.

I need to create a buffer for each blue Area (which I can do) with 5km and then get the amount of people that live in this buffer. Assumption is that people live normally distributed in the green Areas. I already have the amount of people per sqkm.


texts: Raster to CSV with x y cell center coordinates, the original raster cell value, and values from another raster avoiding shapefile

texts:
I am looking for FOSS (maybe GDAL/OGR) solutions to a problem.
I need to go from raster (200 million records) to CSV. It is easy enough to do via point but the point file breaks the allowable 2GB “any component part” shapefile rule so I would technically have to go to another point format such as file geodb but would like to avoid points altogether. I need to get another raster value into the CSV as well.

The CSV must carry over

The original raster Value (16 –bit integer).
An X Y value (technically just a y is needed) that is the center
of the raster cell. This is in DD to double precision.
A value from another raster. Double precision.

What I’ve done.

ArcGIS

Raster (file GeoDB) to Point (file GeoDB).
Add X Y coordiantes to the points.
Extract Value to Points (from the
other raster)
Table (file geodb) to Table (csv).

FOSS

Run ArcGIS to get the raster out of the GeoDB and to tiff.
I cannot find a way to go file geodb raster to tiff without ArcGIS so
as I need to use ArcGIS I just go file geodb raster to point into a file GeoDB point feature class.
OGR 2 OGR (file geodb to csv) export to csv carrying across X
any Y using geometry as WKT.
Still need to get the “other raster
value” into (1) or (2) remember I cannot go to shapefile.
I would like to avoid point altogether.

So I plan to get the file geodatabase rasters into GeoTiff to find a FOSS solution. I would like to go straight from Raster to CSV carrying across the raster value in each cell, generating the x,y value for the raster cell (center), and then get another raster value into the CSV via the x,y CSV coordinate pair. This would avoid going to point and avoid using Arc.


texts: Legend of GeoDataFrame does not sort categorical data values

texts:
I try to plot a map of temperature with geopandas in Google Colaboratory. This is my dataset:

id = '1LOywA8JPP6vOFDP-CBOBgkSmX7zv2pDc'

downloaded = drive.CreateFile({'id':id})
downloaded.GetContentFile('temp_zone.geojson')
print('The files was downloaded')

temp_zone = gpd.read_file('temp_zone.geojson')
print(temp_zone.rango.unique())
temp_zone.head()

And this is the output:

How you can see when I print the column of "rango", this have not any order, so I will sort as I want:

ordered_temp = [ '< 8', '8 - 12', '12 - 16','16 - 20', '20 - 22', '22 - 24','24 - 26','26 - 28', '> 28' ]
temp_zone['rango'] = temp_zone.rango.astype('category',
                                        ordered=True,
                                        categories=ordered_temp)
print(temp_zone.rango.unique())

Finally, when I try to plot the map, the values in the legend are not sorted as wrote before in the categorical data, because the value '<8' is near the last position when this should be in the first position. The image is something like this:

Anyone does know how to fix it?


texts: Fix slivers (Holes) between polygons ( generated after removing spikes)

texts:
I am trying to fix topological errors in this Shapefile of small areas of a state (Bahia) in Brazil. After applying prepair and PPrepair, the resulting shapefile has some spikes:

It looks like a line from far, but when we zoom in we can see that the spike has some area:

Then I removed the spikes with a Postgis buffer function:

SELECT geom,
  ST_Buffer(  ST_Buffer(  ST_Buffer(geom,0.00001, 'join=mitre'),-0.00002,
  'join=mitre'),0.00001, 'join=mitre') AS geom_spike_removed
  FROM repaired_polygons;

resulting in this, (I overlayed the layer) and got some holes between polygons:

How do I remove the slivers so that the empty area is attributed to the nearest polygon?

The spike and subsequent sliver is about 0.071mm in width.


texts:
I'm working on a project which examines different approaches for indoor positioning. All of these approaches return coordinates (lat/lon). What I need to do is create a digital map to compare returned coordinates to real ones, so I can calculate the error/precision of certain methods. How can I do this?

I was thinking, that I could step in the corners of my room and get precise coordinates and use them as comparison, but obviously it is not as easy as it seems. I have seen a lot of bachelor thesis on this topic, but no-one really shows the implementation of digital map. Any ideas how to achieve this?


texts: Getting qgis-ltr-bin.exe - System Error?

texts:
I have installed QGIS 2.8 but when I double click on QGIS desktop it shows an error

  qgis-ltr-bin.exe - System Error

  The program can't start because api-ms-win-crt-runtime-l1-1-0.dll is
  missing from your computer. Try reinstalling the program to fix this
  problem.


texts: Quick way to determine if facing a given lat/lon pair with a heading

texts:
I'm looking for a fast way, preferably without trigonometric operations (just multiply, divide, add, sqrt and subtract) to determine if a point is visible. E.g. say I am facing 0 degrees N, I am at 51.40,-1.08, the camera has a HFOV of 30 degrees and the point is at 51.20,-1.085 it should return TRUE, because the point would be facing me and is within the 30 degree FOV.

It's required that it returns TRUE (false positive) instead of FALSE given edge cases, otherwise points may be hidden when they are actually visible. Currently I have my augmented reality project using atan2 to calculate the actual angle but if I can shave calculations off from points which are definitely not visible every little bit helps. I only need it to be accurate over distances of less than 2 km.

Here's my code to calculate horizontal and vertical position. It's only accurate over <5km distances but accuracy isn't critical anyway as those waypoints aren't drawn. This project is open source, so it's okay for me to post a lot of code.

res->est_dist = hypot2(p_viewer->p.lat - p_subj->lat, p_viewer->p.lon - p_subj->lon) * 2000.0f;
// Save precious cycles if outside of visible world.
if(res->est_dist > cliph)
    goto quick_exit;
// Compute the horizontal and vertical angle to the point.
res->h_angle = angle_dist_deg(RAD2DEG(atan2(p_viewer->p.lon - p_subj->lon, p_viewer->p.lat - p_subj->lat)), p_viewer->yaw);
res->v_angle = RAD2DEG(atan2(p_viewer->p.alt - p_subj->alt, res->est_dist)); // constant scale factor: needs to be adjusted
// Normalize the results to fit in the field of view of the camera.
res->h_angle += p_viewer->hfov / 2;
res->v_angle += p_viewer->vfov / 2;
// Compute projected X and Y position.
res->pos.x = (res->h_angle / p_viewer->hfov) * p_viewer->width;
res->pos.y = (res->v_angle / p_viewer->vfov) * p_viewer->height;

See that call to atan2 for the h_angle? If I could eliminate it by saying "this point is definitely not visible" it would be nice.


texts: What is the best method to smooth TIN Surfaces?

texts:
I have a 3D TIN which shows an avocado sizes in Z-values. The 3D visualization is perfect, however, is not our needs because it requires too much interpretation from our end users perspective.

How can we create continuous surfaces with the same output, while we want to avoid using the interpolation toolbox?

The image below is the TIN surface clipped to a farm boundary.  The lighter green is indicative of size 48 fruits, which are the largest.  Ideally, I would like to create a continuous raster with similar output, such that interpolating the surface will not create an equal distribution across the farm.


texts: Finding movement starting points in QGIS

texts:
I have multiple csv files (24+) with objects categorized by IDs, with their latitude and longitude, and timestamp. Each csv with a different timestamp specified in a column.
There are multiple rows of latitudes and longitudes that often repeat themselves across time, but there are sudden changes in either latitude or longitude and I need to find the specific point where that movement started. I would want to capture even minimal Lat/longs changes, i.e. ID 101 has coordinates consisting of -20.0, 90, and an hour later it is at -20.0001, 90.002.
 I could also combine all the separate csv files into one. Basically, if I were to look at it in table form, I want to find the last row prior to any change in either the latitude or the longitude for each object.
I am not interested in finding the last row for each id if it is not preceded by a new lat/long.  There can be multiple rows with the new lat/longs in the combined file, so remove_duplicates() in pandas via python is not an option.
In other words, if we go the QGIS route for a solution, I would only want the output to have the points prior to a movement in time.

Is there a plugin that would make finding the point prior to coordinate change viable to perform this in QGIS?


texts: Calculating Rumple_Index for clipped Inventory Plot LAS files

texts:
I have successfully implemented the Rumple_Index example provided in the LidR book. I am now trying to calculate the same Rumple_Index for my forest field plots. I am trying to integrate the Rumple_Index function used for the landscape analysis at my plot scale. I am using code to do this that I have used to calculate Zmetrics for the plots already (substituting the rumple_Index function).
The error I am getting now is:
Error in UseMethod("readLAS", files) :
  no applicable method for 'readLAS' applied to an object of class "c('cluster', 'multiprocess', 'future', 'function')"
Called from: readLAS(cluster)

This is a simplified version of the code I am trying to run.
rumple_index_surface = function(cluster, res)
{
  las = readLAS(cluster)
  if (is.empty(las)) return(NULL)
  las <- filter_surfacepoints(las, 1)
  rumple <- grid_metrics(las, ~rumple_index(X,Y,Z), res)
  bbox <- raster::extent(cluster)
  rumple <- raster::crop(rumple, bbox)
  return(rumple)
}

LASlist <- list.files(path to clipped LAZ plots) # create a list of the LAZ files to cycle through

getMetrics <- data.frame()

for ( i in LASlist) {
  lidar = readLAS(i)
  getMetrics <- rbind(getMetrics, cloud_metrics(lidar, ~rumple_index_surface(cluster, res = 25)))
}

How can my code be improved?


texts: List data sources from multiple .mxds contained in one folder using Arcpy

texts:
I have a folder containing multiple .mxd files. I wish to list all of the layer names and data sources for each .mxd from the folder. I am using the following code which lists the above for one .mxd but how do I loop the code to work through all .mxds in the folder?

import arcpy

mxd = arcpy.mapping.MapDocument("D:\MXDs\Air_quality.mxd")
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        print " MIS Name: " + lyr.name + "  Data Source: " + lyr.dataSource


texts: Use gdal to set multiband raster values to no data

texts:
I need to set a 3 band multiband raster (GeoTIff) where band 1 = 0 AND band 2 = 0 AND band 3 = 0 then this is no data (0,0,0 is No Data). Not OR in each band but AND across the 3 bands.

I think it is GDALWARP -dstnodata but am unsure of the syntax.

How do I set the No Data value as 0,0,0 in the header?


texts: Creating instances of QGiS in for loop in standalone python script

texts:
I have a Python script which processes data with geographic information. At some stage of my script I need to use an instance of QGiS to make some shapefiles and do some intersects and so on. For now, I have created an instance of QGiS at the beginning, as it is written to do in documentation, however like this I free the memory and remove data providers and memory layer at the end of my script.

qgs = QgsApplication(sys.argv, True)
qgs.initQgis()

startScripting()

qgs.exec_()
qgs.exitQgis()

It is working, but unfortunately from time to time data from previous processing are also used in next processing.
In startScripting() I have a for loop:

for pro in users[idUser].keys():
    [result] = rep.report(users[idUser][pro],idUser,pro)

And all geographic data are process there (in report()).
So to be honest I need to create instance of QGiS in for loop to have a clear instance to proceed. I have tried something like this:

for pro in users[idUser].keys():
    qgs = QgsApplication(sys.argv, True)
    qgs.initQgis()
    [result] = rep.report(users[idUser][pro],idUser,pro)
    qgs.exec_()
    qgs.exitQgis()

Of course I have deleted creating instance of QGiS at the beginning in this solution.

However if I try creating QGiS instance in for loop, then everything is fine but only for first loop, the second one doesn't start. I think that qgs.exitQgis() stop my whole script.

Maybe I do something the wrong way?
Or maybe there is another way to clear data providers, memory layers, without closing instance of QGiS?


texts: Encoding issue in QGIS while exporting a layer with Cyrillic characters from gpkg to CSV

texts:
I need to export gpkg layer contained Cyrillic characters to a spreadsheet, e.g. csv file. I do that using Export -> Save Features As and pick UTF-8 encoding. I tried encodings: Cyrillic, different UTF's but nothing works. Just UTF 16LE keeps very few Cyrillic words. I also tried to find encoding in Advanced Settings according to an answer from here but it is not available in QGIS 3.16:

Does anyone know how to fix?


texts: How do I get the domain from an OGRSpatialReference?

texts:
I would like to know what the projected bounds are for a spatial reference using OGRSpatialReference.  Is there some way to use the GetProjParm method to get this information?

An example of what I'm looking for can be seen as the "projected bounds" and "wgs84 bounds" on the following spatialreference.org page:

http://www.spatialreference.org/ref/epsg/3786/


texts: Rasterio and OpenCV shows two different pixel arrays for same image

texts:
I am using Rasterio to open a panchromatic GeoTiff image file.
The syntax I used to open the GeoTIFF file using Rasterio:
import rasterio

src = rasterio.open('Birmingham_Shuttleworth_s106_20200409T191236Z_Basic_Panchromatic/20200409_191236_ssc6d3_0015_basic_panchromatic.tif')
band_1 = src.read(1)

The array I can see after doing
In [23]: band_1
Out[23]:
array([[3466, 3865, 3929, ..., 4472, 4520, 4456],
       [3546, 3657, 3721, ..., 4568, 4728, 4504],
       [3562, 3689, 3945, ..., 4696, 4504, 4568],
       ...,
       [5606, 5622, 5478, ..., 4376, 4264, 4424],
       [5558, 5526, 5462, ..., 4600, 4504, 4536],
       [5558, 5510, 5462, ..., 4824, 4696, 5175]], dtype=uint16)

However, when I open the same GeoTIFF file using OpenCV:
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

img = cv.imread('Birmingham_Shuttleworth_s106_20200409T191236Z_Basic_Panchromatic/20200409_191236_ssc6d3_0015_basic_panchromatic.tif',0)

I can see:
In [44]: img
Out[44]:
array([[13, 15, 15, ..., 17, 17, 17],
       [13, 14, 14, ...,  0,  0,  0],
       [13, 14, 15, ...,  0,  0,  0],
       ...,
       [21, 21, 21, ..., 13,  9, 11],
       [21, 21, 21, ...,  0,  0,  0],
       [21, 21, 21, ...,  0,  0,  0]], dtype=uint8)

I see that the pixel type(?) is uint8 in OpenCV and uint16 in Rasterio. Below are some of the questions that I have.

What and why is this happening?
What are those numbers that I see in the arrays? Are they pixel values? If yes, why are they not from 0-255? Do panchromatic image pixels have a different range?

I have just started image processing using Python.


texts: How to draw an overlay using KML on our own app in Android?

texts:
Repost from here

I want to draw an overlay using KML files, just look like below but on my own app in MapActivity.

31.kml, display district in Jakarta, Indonesia
or the example that I get from googling

I've already tried this code, but it just opens Google map app that already gave me a map with KML overlays.

final Intent myIntent = new Intent(android.content.Intent.ACTION_VIEW,
Uri.parse("geo:0,0? q=http://dl.dropbox.com/u/32264286/Unggahan/31.kml"));startActivity(myIntent);

Is it possible and how can I implement this?


texts: Interpolation producing weird result QGIS

texts:
I have some points I want to interpolate.

Output:

Anybody got a good idea on this one?


texts: Getting travel time while performing pgrouting i.e. what are parameters and function required?

texts:
While doing the pgrouting to find the shortest distance I got success in getting the output but when I thought to come up with the time consumed for travelling, I cannot find the proper source to get the time consumed.

If there is a solution, can you let me know about the function to be used for that in pgrouting?


texts:
I am attempting to setup a website on GoDaddy and am having some problems with the OpenLayers API. The API hangs and I end up with a message to continue/debug/stop the script. I've tried both using OpenLayers using a CDN link and copying it to my site, but same results.

To make this question simpler, I went ahead and created a page using the osm.html page which is as simple as one can get. The URL is www.enjoytheparks.com/test.html.

Even this simple page using OpenLayers exhibits the same issue. The specific line is OpenLayers.js:59. I should probably use the non-minified version, but I'm hoping someone has seen this before and it is an easy fix, rather than having to debug through the OpenLayers.js file.


texts: Merging hexagons in layer to create multiple polygons of equal size using QGIS?

texts:
I'm trying to create a hexagon cartogram of electoral constituencies for Pakistan in QGIS. I've got a shapefile of Pakistan's boundaries and I've loaded a layer of overlapping hexagons on top using the MMQGIS plugin. Now I want to merge the hexagons into 272 polygons consisting of equal number of hexagons.

I want to keep the overall shape of the country intact, so any remaining polygons that cannot be added should be lying on the outer boundaries so I can keep the overall shape close to the shape of the country without any gaps in between.

The above example image shows a couple of hexagons manually merged into polygons. I want all 272 polygons to comprise of an equal number of hexagons each. How can I do this in QGIS?


texts: Print composer layout changes when converting to PDF

texts:
I have the problem of the layout/formatting of my QGIS 2.13.2 in the print converter changing whenever I try to print to PDF Form.
In the print converter my layout looks good and the shapefile is at the right place. But as soon as I convert it to PDF Form, my shapefile moves a little bit up.

Is it again a problem of georeferencing?


texts: Ordering layers in OpenLayers and Google styled maps?

texts:
I'm using OpenLayers with Google styled maps. My biggest problem is KML polygon layer which needs to be under roads layer but it isn't (see screenshot).

What am I doing wrong?

This is my coding:

   <title>OpenLayers Google (v3) Layer Example - demo with styled maps, traffic, fusion tables</title>


texts: How to use calculate geometry on non versioned data?

texts:
Every time I attempt to use calculate the geometry of some shapefiles I receive this message: ERROR 130051: Input feature class is not registered as versioned.
I did not understand why a mathematical calculation needed versioning.
Key facts

The data is all stored locally.
The layer that I want measured was created through multiple selections and merges.

The layer started out as two different layers. The first was a map of property parcels in my county. The second was a geological survey for the soil in the county. A selection of large plots was made, then merged with the soil layer so that the individual soil types were divided within each large property.
The goal is to calculate a theoretical crop yield based on soil quality, but there is no way to measure the amount of soil in each parcel because ArcGIS Pro will not make the measurement because the data is not versioned. I cannot publish it as a web layer since selections are forbidden to export. I am not interesting in publishing my unfinished work either. If there is a way to go into the back door to have ArcGIS Pro do this or if I have to use Python please let me know. I do not think this should have such a complex solution but I am all ears.


texts: How to add QGIS Python API Intellisense support in Aptana Studio3/Pydev

texts:
I'm developing QGIS plugins using PyDev in Aptana Studio 3 and I would like to know if it's possible to make it aware of the QGIS API and add some "intellisense" functionnalities when programming. For example, I would like to see the available methods of an object when typing the dot "." after it...

UPDATE:
I also see this message:

  module not in pythonpath

When Hovering any QGIS object in code... Is it supposed to show me something more usefull?

Thanks a lot

Etienne


texts:
I have a PostGIS query which works fine on its own:

select distinct
    gis.filelinks_qtrqtrs_view.entity,
    gis.filelinks_qtrqtrs_view.type,
    gis.filelinks_qtrqtrs_view.url
from
    public.intersection_layer as f,
    gis.filelinks_qtrqtrs_view
where
    st_intersects(f.geom,gis.filelinks_qtrqtrs_view.geom) = true;

I want to create a plpgsql function which will perform this query with the table 'public.intersection_layer' passed to it. The template that I have created:

    CREATE OR REPLACE FUNCTION gis.get_poly_fence(table_name text)
  RETURNS TABLE (entity text,type text,url text)
AS $$
BEGIN
    return query select distinct
        gis.filelinks_qtrqtrs_view.entity,
        gis.filelinks_qtrqtrs_view.type,
        gis.filelinks_qtrqtrs_view.url
    from
        table_name as f,
        gis.filelinks_qtrqtrs_view
    where
        st_intersects(f.geom,gis.filelinks_qtrqtrs_view.geom) = true;
END; $$
LANGUAGE plpgsql;

If I pass a query:

select * from gis.get_poly_fence('public.intersection_layer');

I get an error of "relation "table_name" does not exist. I can't find any way of delimiting the tablename 'public.intersection_layer' argument that works to pass to the 'gis.get_poly_fence' function. I've used functions before (obviously not too many)-- what is going wrong?


texts: Counting points between two intersections: Roadway Network

texts:
I have a roadway network as a line feature and crashes as point features.

I need to count the number of points that lie between two intersections of roadways. The biggest problem is the road is also segmented and is not continuous throughout.

Please look at the below image:

Image 1 shows that I have to count points b/w two intersections.

Image 2 shows roadway segments


texts: Enable groupbox only when checkbox is checked - pyqgis

texts: I have a groupbox which is always enable. I want it to be enable only when user want to edit something, so I made a checkbox for it. However, I have no idea how to 'link' if checkbox if selected with groupbox switching on. I want groupbox to be greyed out when disable.


texts: Is the error message "The Microsoft Jet database engine can not find the input table or query" the result of a hardware limitation?

texts:
I'm using ArcMap 10.2.1.

When I try to compute a union of two shapefiles, I get this message:

  "The Microsoft Jet database engine can not find the input table or query. Make sure it exists and its name is correctly spelled".

The table name that can't be found corresponds to the output shapefile that is given by me (name and location).

I have already performed this operation successfully with two similar files, but now this is the only thing I get.

Can this be a memory problem?


texts: Training Manual 2.8 and exercise data differences

texts:
I am a new user and trying to teach myself about QGIS using the 2.8 Training Manual and the training_manual_exercise_data downloaded from QGIS.

I am using the PDF of the manual, which seems to be consistent with the on-line version.

But even with the very first exercise the shp file it ask you to load have different names or some don't exist. There are no files available called Rivers or Buildings, and the ones there have _335 suffixes to the names. Okay, I can work round this, but the examples are then different.

The next stage is to load landuse.sqlite, which does not appear to exist in the exercise data at all.

So I am wondering if there is a different set of exercise data that is the one I should be using.


texts: Finding boundary co-ordinates from given set of point co-ordinates?

texts:
Given a set of coordinates, How do we find the boundary coordinates.

<== Figure 1
Given the coordinates in the above set, How can I get the coordinates on the red boundary. Boundary is the polygon which is formed by the input coordinates for vertices, in such a way that it maximizes the area.

I am working on an app which searches properties within 'x' miles of a city. What I have is:

Coordinates of all the properties.
A set of coordinates for each city (I have one coordinate for each zip. And since most cities have more than one zip, Every city has a set of coordinates)

The reason I am asking for the maximum area is so that I don't come up with a polygon like the one below:

  <== Figure 2

What I need is the algorithm to come up with the set of coordinates for the boundary. An algorithm which will allow me to come up with boundary coordinates for Figure 1.


texts:
I am using the OSM open Roads WMS layer that comes included with MapInfo. I have noticed that when I load this layer the projection of the Map window changes to 'Popular Visualization CRS.'

All the other layers I use are set to British National Grid. Is there a way to change the projection of the WMS layer so it matches British National Grid?


texts: QGIS Layout: Making a map frame in a certain scale bar style with world coordinates of the map (view) extent

texts:
Currently I am making some layouts for my some templates we are going to use. We would also like to have a certain border "frame" on to make the coordinate grid clear. The frame is supposed to look like this:

At first this looks simple we just need to make the frame based on a certain scale bar style, but here things get a bit more complex only that but they also have integrated the X and Y Coordinates of the map (view) extent. Sometimes like this for certain maps:

But other times like this:

So how do we do this? Not to mention this has to update dynamicly since we are making "templates" layouts, this has to automaticity update and soon as the view and extent changes.


texts: How to determine column value based on lowest value from other column?

texts:
I have the following dataframe with same value in 'nearest_beacon' column but different distance in 'vms_distance' column.

My Dataframe:

nearest_beacon  vms_distance associated
4548780      0.486456        vms
4548780      0.468065        vms
4548780      0.337609        vms
4548780      0.363601        vms
4548780      0.210566        vms
4548780      0.197327        vms*
4548780      0.285390        vms
4548780      0.216443        vms
4548780      0.441454        vms
4548780      0.337533        vms

I want to determine the 'associated' column just one row (*) which has 'vms' value considering it has the lower value in 'vms_distance' column, and the rest is 'no_vms'.

Expected Result:

nearest_beacon  vms_distance associated
4548780      0.486456        no_vms
4548780      0.468065        no_vms
4548780      0.337609        no_vms
4548780      0.363601        no_vms
4548780      0.210566        no_vms
4548780      0.197327        vms
4548780      0.285390        no_vms
4548780      0.216443        no_vms
4548780      0.441454        no_vms
4548780      0.337533        no_vms


texts: "Spatial join" between raster and shapefille- give raster the shapefile values with python

texts:
I have raster which has NDVI values .In addition, I have shapefile(geopandas) on top of this raster but it doesn't fully covers the raster (illustration):

the shapefile has column which contain data about each polygon:

I would like to  do something like spatial join between them, to have this data (e.g X,Y,Z,A,C) stored for each pixel,
so in the end i'lll have "new band" which has one band of NDVI and one with this letters data,
so for example, pixel that doesn't fall inside the polygon gets value "null" or "nodata", and pixels that falls inside the polygon will get the value of the letter according to where they fall inside (if thery are inside polygon X they get value X, if they are inside polygon Z then they get Z ect.
yet I haven't found any way to do that, I though maybe to clip the raster with the shapefile and then give it the values but also then i'm not sure how to give it the different letters values.
My end goal: To have the shapefile data (which is stores in one column in the table) stored for each pixel, if the pixel is not inside the polygon - give no data value. Kind of spatial join for raster-pixel .
Edit: i'm using python with jupyter notebook .
I had the thought maybe to extract the raster to point and then try to do spatial join and back to raster but it's huge number of pixels
Edit2:
I have tried Jose's solution to use this function:
def rasterise_me(raster, vector, attribute,
                fname_out="", format="MEM"):
    """Rasterises a vector dataset by attribute to match a given
    raster dataset. This functions allows for the raster and vector
    to have different projections, and will ensure that the output
    is consistent with the input raster.

    By default, it returns a handle to an open GDAL dataset that you
    can e.g. `ReadAsArray`. If you want to generate a  GTiff on disk,
    set format to `GTiff` and `fname_out` to a sensible filename.

    Parameters
    ----------
    raster: str
        The raster filaname used as input. It will not be overwritten.
    vector: str
        The vector filename
    attribute: str
        The attribute that you want to rasterize. Ideally, this is
        numeric.
    fname_out: str, optional
        The output filename.
    format: str, optional
        The output file format, such as GTiff, or whatever else GDAL
        understands
    """
    # Open input raster file. Need to do this to figure out
    # extent, projection & resolution.
    g = gdal.Open(raster)
    geoT = g.GetGeoTransform()
    nx, ny = g.RasterXSize, g.RasterYSize
    srs = g.GetProjection()
    min_x = min(geoT[0], geoT[0]+nx*geoT[1])
    max_x = max(geoT[0], geoT[0]+nx*geoT[1])
    min_y = min(geoT[3], geoT[3] + geoT[-1]*ny)
    max_y = max(geoT[3], geoT[3] + geoT[-1]*ny)
    # Reproject vector to match raster file
    vector_tmp = gdal.VectorTranslate("", vector, format="Memory",
                                    dstSRS=srs)

but when I try to apply this I get error:
    # Do the magic
    ds_dst= gdal.Rasterize(img, plots, attribute='Crop',
                        outputSRS=srs, xRes=geoT[1], yRes=geoT[-1],
                        outputBounds=[min_x, min_y, max_x, max_y],
                        format=format, outputType=gdal.GDT_Int32)
    return ds_dst

NameError: name 'srs' is not defined


texts:
I have QGIS 3.16.3-Hannover installed with Qt 5.11.2 and 64-bit Windows 10. In two more computers I also have it installed in similar conditions.
I am using the PyQGIS to stand-alone applications from QGIS and I think there is something that is not working well within PyQGIS because I cannot use the QWebEngineView and the QWebKit yes.
The symptoms I have are:
1.- If in the QtDesigner that is installed with QGIS, I add the QWebEngineView widget in the design, the program directly crash and stops running.
2.- If in a stand-alone application I start a QWebEngineView it gives me an error:
Code:
from qgis.PyQt.QtWidgets import QWidget, QApplication, QMainWindow, QFileDialog
from PyQt5.QtWebEngineWidgets import QWebEngineView

import sys

class temp (QWebEngineView):
    def __init __ (self):
        super () .__ init __ ()

if __name__ == '__main__':
    # We create the instance of the created class and show the window on the screen
    app = QApplication (sys.argv)
    win = temp ()
    win.show ()
    sys.exit (app.exec_ ())

and the error it gives is:
Qt WebEngine ICU data not found at C: /OSGeo4W64/apps/Qt5/resources.Trying parent directory ...
Qt WebEngine ICU data not found at C: / OSGeo4W64 / apps / Qt5. Trying application directory ...
Qt WebEngine ICU data not found at C: / OSGeo4W64 / apps / Python37. Trying fallback directory ... The application MAY NOT work. Installed
Qt WebEngine locales directory not found at location C: / OSGeo4W64 / apps / Qt5 / translationss \ qtwebengine_locales. Trying application directory ...
Qt WebEngine locales directory not found at location C: / OSGeo4W64 / apps / Python37 \ qtwebengine_locales. Trying fallback directory ... Translations MAY NOT not be correct.
[0205 / 214808.373: ERROR: icu_util.cc (179)] Invalid file descriptor to
ICU data received. [0205 / 214808.491: FATAL: icu_util.cc (311)] Check failed: result. Backtrace:
``
I have looked at this link:

https://issues.qgis.org/issues/18155
3 years ago it seems that this error could not be solved ...

Then with this other:

https://forum.qt.io/topic/70457/qt-5-7-qtwebenginewidgets-not-working-on-linux/3

It seems that the following folders are missing, among others, that I do not have in my installation: resources, qtwebengine_locales, ...

I have uninstalled the qwebengine with the OsGeo4W setup. And then I have reinstalled it.

I have seen that QGis uses version 5.10.1, I have tried to install it through pip install, but it tells me that there is only the version from 5.12.

I think it must be that the QWebEngineView is not installed correctly with QGIS. Because if that works if I use it with the installation of python38 or python39.


texts: How can I improve multi-core processor utilization % in Mapinfo?

texts:
What would be the best way to get large tasks to run in MapInfo without booking all afternoon for a few split and merge queries? I have two quad-core processors in a Win7 machine, but mapinfo will only use a single core, 13% of my total processing power.

Is there any way to get MapInfo use multiple cores, or just to use more of my processing power?
If not, do other packages like QGIS or ESRI work better?

I work with very large polygons where most tasks are scripted and scheduled overnight. When I'm asked a non-standard question (Like right now) I spend an hour+ on every step. Is this just my penance for having a very simple file structure with very large files?

EDIT: In the future, I would make a few copies of my input files and try running multiple parts of the query on the  different sets of files in different MapInfo instances, as per Peter's suggestion below. I would also probably use 'Enclose' to get the extents of my polygon files rather than 'Joining' all the often overlapping files together to get information about area, and at the same time run a parallel MapInfo instance for separate queries about population and then join the output at the end stage.
However, this is only a partial solution, as I would still like faster processing in order to be able to try and fail quickly, and to get an output in a reasonable time no matter how silly and inefficient my workflow is. And a pony too, I suppose.


texts: How to import .csv file to CartoDB using SQL?


texts: Is it feasible to reclassify a raster which has been resampled/smoothed for display purposes?

texts:
Background:

I have a Google Earth KMZ file which shows ground movements after an earthquake, arranged into discrete classes:

From the KMZ file I'm able to download the original PNG images, which are multi-band (RGBA). I used the mosaic tool to create a single, seamless raster over my area of interest. I need to reclassify this image to match the legend above.

The problem:

The problem is that while there are 13 discrete classes in the legend, there is a far greater range of values in the mosaicked raster.

I believe this has occurred due to resampling and smoothing of the images (which occurred before the data provider gave them to me) which has resulted in interpolation between the known values. Here is a screenshot of the cells at close magnification:

For example, the yellow in the legend corresponds to RGB values of 255,255,0 but the yellows in the mosaicked raster vary considerably (eg 239,239,15 or 246,246,8, etc).

My question:

Is it feasible to "resurrect" the mosaicked image, and convert it into something approximating the values in the legend?

I realise that this approach isn't ideal, and that it would be far better to start with the original source dataset - but this isn't possible in this case.


texts: Why latitude/longitude and not the other way round?

texts:
When you google "latitude/longitude" you get 10 times more results than "longitude/latitude".  This I find very confusing as "x/y" seems 20 times more common than "y/x".  But on a map Latitude is on the Y axis and longitude on the X.

Maybe I'm just venting, but maybe I'm missing a clue that makes all this madness logical.

Any explanations?


texts: Getting the raster file of elevation UK with 1 km grid?

texts: I am trying to get raster file of elevation UK with 1 km grid? The EU-DEM has an elevation of 30 meters and 25-meter grids. I could not find any source for 1 km.


texts: QGIS One Shapefile, Two attribute colums, both need to be styled in dependence of each other

texts:
I didn´t find an easy solution for my problem.
Im working with QGIS and i got one shapefile but i want to implement two atributes together.

I need to make a bubble with different sizes based on attribute column 1 (Values 0-2) and
the color in the bubble Should to be adjusted based on attribute column 2 (Values from 0 - 50, 5 colors (0-10=red, 10-20=yellow, etc.).

The color in the bubble should fill out the whole bubble. Can you help me?


texts: Name of a custom Goode homolosine projection

texts: What is the specific name of this homolosine projection? I see that America is closer to Europe/Africa than in the typical homolosine.


texts:
I want to create a windmap. So I created a csv with coordinates with values for wind speed and direction, these represent my weather stations. I loaded the csv up in QGIS and marked a station. Then I chose Raster ==> interpolation ==> I chose the value direction and clicked "add". Then I chose the interpolation method Inverse Distance Weighting (IDW) for processing. Now I have a raster for my speed value and I have also processed one for the direction value.

Here is my Question:

Now I want to put a regular points grid over my interpolated rasters but it doesn't shows up? So I go to Vector ==> research tools ==> regular points.

Every time I want to create a regular points grid it proceeds and gives a result wich is an empty "regular points" layer?

I tried to use a layer, as a mask, and also tried it without a mask layer.

There are no errors or obvious failures? I tried it with other rasters and there were no problems at all.

It would be great to get some ideas and help.


texts:
I'm using OL4 on my project and I have a feature layer where I can click and a popup appears, I also use a draw layer where I can plot a Marker (in the example it's just a point)

My problem is that if I plot the marker above the feature layer, the click propagate through the draw layer till the feature layer and the popup appears, I know it exists the stopPropagation() event but I really don't find the correct place to insert it.

I tried to insert the event on marklayer and map on(click...) but doesn't work.

Here my code:

var container = document.getElementById('popup');
var content = document.getElementById('popup-content');
var closer = document.getElementById('popup-closer');
var overlay = new ol.Overlay({
    element: container,
    autoPan: true,
    autoPanAnimation: {
      duration: 250
    }
});
    closer.onclick = function() {
      overlay.setPosition(undefined);
      closer.blur();
      return false;
    };

// DRAW Layer
var Msource = new ol.source.Vector();
var markLayer = new ol.layer.Vector({
   source: Msource,
   style: new ol.style.Style({
     fill: new ol.style.Fill({
     color: 'rgba(255, 255, 255, 0.2)'
   }),
   stroke: new ol.style.Stroke({
     color: '#328cc1',
     width: 4
   }),
   image: new ol.style.Circle({
     radius: 7,
     fill: new ol.style.Fill({
       color: '#328cc1'
     })
   })
   })
 });

// DRAW Layer
var Msource = new ol.source.Vector();
var markLayer = new ol.layer.Vector({
   source: Msource,
   style: new ol.style.Style({
     fill: new ol.style.Fill({
     color: 'rgba(255, 255, 255, 0.2)'
   }),
   stroke: new ol.style.Stroke({
     color: '#328cc1',
     width: 4
   }),
   image: new ol.style.Circle({
     radius: 7,
     fill: new ol.style.Fill({
       color: '#328cc1'
     })
   })
   })
 });

// MAP
      var map = new ol.Map({
        layers: [
          new ol.layer.Tile({
            source: new ol.source.TileJSON({
              url: 'https://api.tiles.mapbox.com/v3/mapbox.natural-earth-hypso-bathy.json?secure',
              crossOrigin: 'anonymous'
            })
          })
        ],
        overlays: [overlay,markLayer],
        target: 'map',
        view: new ol.View({
          center: [0, 0],
          zoom: 2
        })
      });

// POPUP INTERACTION
map.on('singleclick', function(evt) {
  var coordinate = evt.coordinate;
  var hdms = ol.coordinate.toStringHDMS(ol.proj.transform(
          coordinate, 'EPSG:3857', 'EPSG:4326'));

  content.innerHTML = '
You clicked here:
' + hdms + '';
  overlay.setPosition(coordinate);
});

// MARKER INTERACTION
$("#marker").click(function(e) {
  addMark("Point");
});

var mark;

function addMark(Type) {
  mark = new ol.interaction.Draw({
    source: Msource,
    type: Type
  });

    map.addInteraction(mark);

    mark.on("drawend", function(){
      //to do stuff
    });

    markLayer.on("change", function(){
      // remove the interaction after you have plotted a marker
      map.removeInteraction(mark);
    });
}

here my CodePen: https://codepen.io/sebalaini/pen/gvPLBJ

On the live system, I have a map layer, feature layer and a drawing layer so not all the map is clickable but, for the purpose of showing my problem, the example is OK, you can see that if you draw a circle you open also the popup


texts: Issue with creating qgs project file with delimited text data

texts:
I have a certain workflow that I would like to automate

I import a geotiff file into QGIS workspace.
I import a csv file into QGIS workspace with contains some x,y,w and h data, already converted to crs of geotiff image.
After that I change the marker type to ellipse, set symbol width and height to w and h(Units -> Map units), set some color and layer opacity.

I want to automate this process using PyQGIS. So far I have managed to create QGIS project and import the geotiff file. But the delimited text file is always failing to load and PyQGIS is reporting the layer as invalid.
import os
import sys
from urllib.parse import unquote, urlencode
from qgis.core import QgsProject, QgsRasterLayer, QgsApplication, QgsVectorLayer, QgsCoordinateReferenceSystem
import gdal,osr

strProjectName = "my_project.qgs"
qgs = QgsApplication([], False)
qgs.initQgis()

project = QgsProject.instance()

ds=gdal.Open("RELATIVE_PATH_TO_TIF")
proj = osr.SpatialReference(wkt=ds.GetProjection())
prj = "EPSG:"+proj.GetAttrValue('AUTHORITY',1)
crs = QgsCoordinateReferenceSystem(ds.GetProjection())
layer = QgsRasterLayer("RELATIVE PATH TO TIF", "ImageLayer")
project.instance().setCrs(crs)
project.addMapLayer(layer)
layer.setCrs(QgsCoordinateReferenceSystem(int(proj.GetAttrValue('AUTHORITY',1)), QgsCoordinateReferenceSystem.EpsgCrsId))

str_csv = "RELATIVE_PATH_TO_CSV?encoding=%s&crs=%s&delimiter=%s&xField=%s&yField=%s" % ('UTF-8',prj.lower(),",","map_x", "map_y")
str_csv = r"{}".format(str_csv)
print(str_csv)
layer_csv = QgsVectorLayer(str_csv,
                           "old_record","delimitedtext")
layer_csv.setCrs(QgsCoordinateReferenceSystem(int(proj.GetAttrValue('AUTHORITY',1)), QgsCoordinateReferenceSystem.EpsgCrsId))
if layer_csv.isValid():
    project.addMapLayer(layer_csv)
    names = [layer.name() for layer in project.instance().mapLayers().values()]
    root = project.instance().layerTreeRoot()
    for name in names:
        layerId = project.instance().mapLayersByName(name)[0]
        myLayerNode = root.findLayer(layerId.id())
        myLayerNode.setExpanded(True)
    project.write("./"+strProjectName)
else:
    print("invalid layer")

I am always getting invalid layer error.
I am running PyQGIS via Anaconda. Python version is 3.6.11 and QGIS version is 3.10
gdal version is 3.0.4.


texts: QGIS Processing algorithm cannot deal with special characters in input

texts:
My aim is to write an algorithm, where the user selects a settlement name, and numeric value, from which a query is generated, and a selection proceeds on the selected layer, than the selection result is exported to a new shape. It works, but when the input settlement name includes a special character, it fails.

I have the following code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

##[Example scripts]=group
##lyr=vector
##nev=string Szeged
##depth=number 3000

from qgis.core import *
import processing

input = processing.getObject(lyr)

print depth
print name
l = input
expr = QgsExpression(" \"MELYSEG\">'{}' AND \"TLP\"LIKE'{}'".format(depth, nev))
it = input.getFeatures( QgsFeatureRequest( expr ) )
ids = [i.id() for i in it]
input.setSelectedFeatures( ids )

QgsVectorFileWriter.writeAsVectorFormat(input, "D:/Data/output.shp", "utf-8", input.crs(), "ESRI Shapefile", 1)

The field type of "MELYSEG" is a float the "TLP" is a string.

When I have an input like this (note the ő character in the name):

I get the following error:

'ascii' codec can't encode character u'\u0151' in position 4: ordinal not in range(128) See log for more details

Any ideas to solve this?


texts:
I am developing a WebGIS application with:

Openlayers 2.12
PostgreSQL 9.2 with PostGIS 2.0 extension
Geoserver 2.2.3
Apache 2.0 configured with Tomcat 7.0.35

PROBLEM

One of the services of the Web application is the generation of an isochronal curve requiring geoprocessing. Here is more detail on the description of the geoprocessing :

Input: coordinated X,Y in " Spherical Mercator " ( 900913 ) of a
point selected by the user on the map and an image (road network
raster).
The geoprocessing : the geoprocessing is applied to the road network
raster and sends back a new image output. This geoprocessing is
already implemented on my machine with Python (without Apache nor
geoserver connections). I use the GDAL library to convert my images
to array onwhich I apply my geoprocessing.
Output: an image which must be displayed on my Web application (a '
layer '  )

SOLUTIONS ?

Do you have a solution to make this service? Which tool or extension should I use?
Here is what I found on the Internet :

The idea would be to integrate my python script (developed on my machine alone) on geoserver. This script should be able to use an image of the road network (hosted on the geoserver) and location (X,Y) of a point sent by the user. Then the script has to answer an image and display customer side with OpenLayers.

WPS geoserver :
this one uses functions predefined at the geoserver level while I want to create my own function !!!

Visual render geoserver side :
the problem is that it does not make call to a script python (but rather in SLD transformation).

Does the python extension of the geoserver allows to generate this kind of geoprocessing? If yes, how to make concretely? I downloaded it on this page http://gridlock.opengeo.org/geoserver/2.2.x/community-latest/  (geoserver-2.2-SNAPSHOT-python-plugin.zi p).
How to exploit(run) the Geoscript library and the Script Hook?

The other solution would be to configure python with Apache. To launch the script python at the Apache level, the result of my script would feed a table in postgres (raster) which is in WMS side geoserver.

EDIT

HERE AN ILLUSTRATION TO UNDERSTAND BETTER THE PROCESS


texts: How can you conditionally halt an FME translation?

texts:
I have an FME workspace that takes a set of street centerlines and converts it into address points. I have been working my way across the state of Missouri with this translator one county at a time, since each county can take up to an hour to process. I select a county at the start of the translator as a choice parameter.

My problem is that I keep having to double check to make sure of what counties I have processed so I do not process a county twice. When I have a file of several million address points, I do not want to be removing duplicates.

What I would like to have FME do is take the county I have selected, and check against a shapefile of counties to see if the county has been processed yet. I have a flag attribute, "COMPLETE" that is set to 0 if the county has not been processed and 1 if it has.

The county shapefile has nothing to do with any other readers or writers in the translator.

Edit: To clarify, the counties are not being processed at all by the translator, only the centerlines.

Since there are over 2 million centerlines, I do not want to read the centerlines at all (much less maintain an attribute on the centerlines) if the county has already been processed.

Edit:
Here is the solution. I moved up the st_county10 reader (the one in the image) to the top of the reader list. This executes before anything else. By using the feature holder, I avoid writing over the st_county10 shapefile if the translation is terminated.

One last piece, since I had a python shutdown script, I had to add this line:

import pyfme
if FME_Status == True:

To the front of the shutdown script to stop execution of the shutdown script too.


texts: Why is using two imbricate UpdateCursor not working?

texts:
I hope my question is not redundant.

I wrote a script to Update a table using two imbricate cursors. The goal was:

to read a records, update some of its values
Make a selection based on this record (e.g. closest neighbors)
Update the neighbors

I have first tried using two imbricate update cursors like the following:

rows = gp.UpdateCursor(Infile)
row = rows.Next()
while row:
    # Read InFile attribute
    val = row.GetValue(Attribute)
    #  Update some attributes
    val = row.SetValue(Field, function(val))
    # Select Neighbors
    gp. gp.GANeighborhoodSelection(Infile, Layer, Condition related to val)
    //.....//
    # Update Neighbors
    urows = gp.UpdateCursor(Infile)
    urow = urows.Next()
    while urow:
        #  Update some attributes
        val = urow.SetValue(Field, function(val))
        urows.UpdateRow(urow)
        urow = urows.Next()
    del urow, urows
    //.....//
    row.UpdateRow(row)
    row = rows.Next()
    //.....//

Using such a code, the second update cursor was calculating correctly (i printed some intermediate results) but the calculation was not applied to the  table (ie. not physical updation).

I modified the code by replacing the second update cursor by a Select by attributes followed by a calculate field. Again, the calculations were not applied to the table. I think there is a lock problem somewhere but I did not found how to solve the problem.

Finally, I replaced my first UpdateCursor by a SearchCursor and then used an UpdateCursor based on a condition in order to reduce the number of record processed. The code look like this:

rows = gp.SearchCursor(Infile)
row = rows.Next()
while row:
    # Read InFile attribute
    val = row.GetValue(Attribute)
    #  Update some attributes
    val = row.SetValue(Field, function(val))
    # Select Neighbors
    gp. gp.GANeighborhoodSelection(Infile, Layer, Condition related to val)
    //.....//
    # Update Neighbors
    urows = gp.UpdateCursor(Infile, Condition based on val)
    urow = urows.Next()
    while urow:
        #  Update some attributes
        val = urow.SetValue(Field, function(val))
        urows.UpdateRow(urow)
        urow = urows.Next()
    del urow, urows
    //.....//
    row.UpdateRow(row)
    row = rows.Next()
    //.....//

Using such a structure, all is working fine. However, calculation is longer (around 5 sec. for 1000 records).

Do you have any idea:

why using 2 imbricate UpdateCursor is not working; and
how to use select by attributes and calculate field within a cursor to fasten the calculation?

I'm using Python 2.5 and ArcGIS Desktop 9.3


texts: How to download hospitals and restaurants from OSM in shp format?

texts:
I need information on hospitals and restaurants in Guatemala.

How can I download this information from OpenStreetMap?

I want to use the data in Arcgis 10.1


texts:
What is the maximum number of records returned by WFS protocol? What is maximum number of records returned by a Feature Layer in ArcGIS Server 10.1? What is maximum number of records the ArcGIS JavaScript 3.1 feature layer can display?

I am facing a peculiar problem in ArcGIS Server 10.1 and ArcGIS JavaScript API 3.1.
In the capability tab, if I set the maximum number of records to more than 13,000 records, then the layer is not displayed. I am using a relational database and joining 5-6 tables to get the desired data.


texts: Print XYZ layer with MapFish Print plugin

texts:
I have a XYZ layer(i hope i call it right) in GeoExt/gxp application:
,{source: "ol",group: "background",fixed: true,type: "OpenLayers.Layer.XYZ",args:
["ЕЭКО","http://maps.rosreestr.ru/ArcGIS/rest/services/BaseMaps/BaseMapAnno/MapServer/tile/$\{z\}/$\{y\}/$\{x\}",{visibility: true}]}

In MapFish documentation says that print plugin can print TMS and Tile Cache. But its print only WMS layers. I gonna set url to XYZ layer or something?
UPDATE
I found working XYZ print in here http://map.geo.admin.ch/. But There is something strange with layer which set to print.
UPDATE2
More interesting that GeoExt.PrintProvider not send XYZ layer to print plugin. But if i give to it gwc layer i can see gwc layer to post parameters:

But this layer not print.
Why PrintProvider not send XYZ layer in post parameters?


texts: Go through items in a composer legend (QGIS 2.6)

texts:
I'm trying to develop a plugin in QGIS 2.6, which opens a composer, add a map and a legend, and (among other things) adjusts the legend.

Is there a way to go through items in a composer legend with PyQgis?
What I'd want in the end is to change the item name, depending on the type of item (group, raster ...).

Here is my code:

for legend in self.composition.items():
    if isinstance(i,QgsComposerLegend):
        legend.setAutoUpdateModel(True)
        legend.setLegendFilterByMapEnabled(True)
        for item in legend.items():
            if item.itemtype() == QgsComposerLegendItem.GroupItem:
                groupItem.setUserText() = 'another text'

Qgis answers that QgsComposerLegend has not attribute called "items" (which I expected).
I looked for a property to iterate on legend items, and I have not found it in the API documentation nor in web searches.


texts: Linking data from Microsoft Access to ArcMap 10.1

texts:
I'm trying to link data in a Microsoft Access file to my work in ArcMap 10.1.
Instructions online show people using the 'Customize' menu option and then selecting 'ArcCatalog' and from there selecting an add 'OLE' function; however, in the list of customize options I do not have an 'ArcCatalog' option to select.

Why might this be and what might I do in order to be able to link the Microsoft access data?

The data set is the HWSD world soils data base.


texts: Calculating last record in feature class using ArcPy?

texts:
I want to automate a task in ArcGIS using arcpy.  The task starts with creating points with 50 meter intervals along a line. Then use a loop to iterate and calculate values 50 meters from 0 to last point.  The problem is I want to calculate the actual distance value of last point and I can't get last record. I try to use rowpt[-1] and use this sample but i don't know how to update the value.  How can I get the last record?

import arcpy
arcpy.env.workspace = r"F:\test\env.gdb"
arcpy.env.overwriteOutput = True
line_lyr = "line"
pt_lyr =  "my_point"
interval = 50
arcpy.CreateFeatureclass_management(r"F:\test\env.gdb",pt_lyr,"POINT")
arcpy.AddField_management(pt_lyr,"Interval","SHORT")
insertCursor = arcpy.da.InsertCursor(pt_lyr, ["SHAPE@XY"]) # this is the pre-existing pt feature class
with arcpy.da.SearchCursor(line_lyr, ['OID@','SHAPE@']) as searchCursor: # this is the line feature on which the points will be based
    for row in searchCursor:
        lengthLine = row[1].length # grab the length of the line feature, i'm using round() here to avoid weird rounding errors that prevent the numberOfPositions from being determined
        if lengthLine % interval == 0:
            numberOfPositions = int(lengthLine // interval) - 1
        else:
            numberOfPositions = int(lengthLine // interval)
            if numberOfPositions >= 0:
                distancezero = 0
                xPoint = row[1].positionAlongLine(distancezero).firstPoint.X
                yPoint = row[1].positionAlongLine(distancezero).firstPoint.Y
                xy = (xPoint, yPoint)
                insertCursor.insertRow([xy])
                for i in range(numberOfPositions):
                    distance = (i + 1) * interval
                    xPoint = row[1].positionAlongLine(distance).firstPoint.X
                    yPoint = row[1].positionAlongLine(distance).firstPoint.Y
                    xy = (xPoint, yPoint)
                    insertCursor.insertRow([xy])
                distanceend = lengthLine
                xPoint = row[1].positionAlongLine(lengthLine).firstPoint.X
                yPoint = row[1].positionAlongLine(lengthLine).firstPoint.Y
                xy = (xPoint, yPoint)
                insertCursor.insertRow([xy])
            del insertCursor
with arcpy.da.UpdateCursor(pt_lyr,"Interval") as ptcursor:
    number = 0
    for rowpt in ptcursor:
        rowpt[0] = number
        number = number + 50
        ptcursor.updateRow(rowpt)
    if rowpt[-1]== number-50:
         rowpt[-1] = 100
         ....?


texts:
QGIS 2.12.1 under Linux Mint 17.2 Cinnamon: I have no touch zoom and pan icon.  Touch screen works for panning and other stuff but I cannot pinch to zoom.

Any fix?


texts: Query all records from ArcGIS Server REST table

texts:
We have a data table stored as table on a ArcGIS server. The data is accessible via REST calls.

I would like to make queries that return all columns from the table, returned as JSON.

https://test.test.org/arcgis/rest/services/testService/MapServer/10/query?where=OBJECTID+%3E0+AND+OBJECTID+%3C+1000

Using the query above I am requesting the first 1000 records from the service. The problem is I don't get all the columns in the table. How could I structure my query to also return all 20+ fields in the table?


texts: How do you use the parent algorithm tools in Extract by location process modeler algorithm?

texts:
I am trying to get the Extract by Location algorithm to work in the process modeler but I can not figure out which algorithms to include in the parent algorithm in order to get it to work.

I've been able to get it to work in the past but going through every combination with 6 different algorithms (64?) is extremely time consuming and isn't working.

The purpose of the model is to extract changes in road network by comparing historical data with new data.  It builds a buffer around the new data and checks to see if the historical data falls outside of the buffer.  Any changes are extracted to be used for field research.  The interstate buffer is used to remove any motorway_links that are not on or near an interstate.


texts: ArcGIS JavaScript API Buffer Units from Dropdown

texts:
I am trying to create a buffer that will have the distance and units inputted by the user.

I'm having trouble getting the JavaScript to recognize the which Unit of Measurement the buffer is being set to. When I click on the map to generate the buffer nothing shows up. This is the JavaScript I wrote:

function doBuffer(evt) {
            map.graphics.clear();
            var params = new esri.tasks.BufferParameters();
            params.geometries = [ evt.mapPoint ];

          //buffer in linear units such as meters, km, miles etc.
            params.distances = [ document.getElementById("distance").value];
            params.unit = document.getElementById("units").value;
            params.outSpatialReference = map.spatialReference;

            gsvc.buffer(params, showBuffer);
        }

This is the HTML I wrote:

Enter the unit and distance to set the buffer.

Distance

                <input id="distance" type="number" name="Distance">

Unit of Measurement

                <select id="units">
                    <option value="esri.tasks.GeometryService.UNIT_FOOT">Feet (ft)</option>
                    <option value="esri.tasks.GeometryService.UNIT_METER">Meters (m)</option>
                    <option value="esri.tasks.GeometryService.UNIT_KILOMETER">Kilometers (km)</option>
                    <option value="esri.tasks.GeometryService.UNIT_STATUTE_MILE">Miles (mi)</option>
                    <option value="esri.tasks.GeometryService.UNIT_NAUTICAL_MILE">Nautical Miles (NM)</option>
                </select>

How can I get my code to recognize which unit of measurement I am setting the buffer to so it can load correctly?


texts: Leaflet photo plugins library compatiblity issue

texts:
I am using the Leaflet photos plugin, it runs perfectly using the Leaflet v0.7.7
but using the newest version of Leaflet v1.4.0
it doesn't work and it gives the error

Uncaught TypeError: Cannot read property 'trim' of undefined
    at h (leaflet.js:5)
    at u (leaflet.js:5)
    at e.on (leaflet.js:5)
    at e.initialize (leaflet.markercluster.js:6)
    at e.initialize (Leaflet.Photo.js:63)
    at new e (leaflet.js:5)
    at Function.L.photo.cluster (Leaflet.Photo.js:80)
    at index.html:65

the web map code is the following

var map = L.map('map').setView([35.924629,5.301136], 14);

var Base =L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: '&copy; DAR TOTAL SOLUTIONS contributors'
}).addTo(map);
var photos = [{
"lat":35.927689,"lng": 5.31083,"name":"","url":"jpg/20170227_084620.jpg","thumbnail":"thumbnail/20170227_084620.jpg"},{
"lat":35.927689,"lng": 5.31083,"name":"","url":"jpg/20170227_084637.jpg","thumbnail":"thumbnail/20170227_084637.jpg"},{
"lat":35.927689,"lng": 5.31083,"name":"","url":"jpg/20170227_084639.jpg","thumbnail":"thumbnail/20170227_084639.jpg"}];
var photoLayer = L.photo.cluster().on('click', function(evt) {
  var photo = evt.layer.photo;
  var template = '
{name}
';
  evt.layer.bindPopup(L.Util.template(template, photo), {
    className: 'leaflet-popup-photo',
    minWidth: 800
  }).openPopup();
});
photoLayer.add(photos).addTo(map);
L.control.layers({
'OSM' :Base,
},{
'Photos':photoLayer,
}).addTo(map);


texts: What Application for Vector and Cartesian Coordinate System

texts:
Goal:
Visualize a picture in SQL server 2012 spatial data based on geometry

Problem:
In order to display a vector picture I need to retrieve X, Y coordinate (cartesian coordinate system) and display it in the SQL server spatial data.

I'm looking for a tool that can draw vector graphic. Every object of line, point or zone you make, you shall enable to retrieve X and Y coordination from the start to the end of the object. For instance line object has start point (0,0) and end point (5,3).

The reason why I need this X,Y data is to use them in SQL server to visualize a vector picture.

  What tool do I need to use?

Other Important Info:

Please take account that I do not want to use geography data because the x-and y-coordinate data that are not defined in relation to any
particular model of the Earth, it's using the geometry datatype with
SRID = 0
The X and Y data retrieve from the special application shall be enable to use in SQL server. In other word, same cartesian coordinate system without measurement of unit.
It's important that this application can apply gif or jpeg picture to act as a background and om top on it, you draw a vector geometry data.
This application should also take account to Well-known text to be customized with SQL server and its spatial data. In other word, take account to OGC standard
The picture from a paper map or jpeg will be a subway map that is a linear cartogram


texts: Calculate sum of points within buffer (ArcMap)

texts: I've created buffers around home origins and am looking to calculate the sum of shapefiles (in this case crime) within each buffer. How would I go about doing this? I want to know how many crimes are within each buffer. Thanks!


texts: Determining what lines have changed in updated shapefile using QGIS

texts:
I am using my county's GIS data to update the roads in OpenStreetMap. I have a copy they published last year and one that was just published, and I would like to find all the LineStrings that have either had their attributes or their geometry changed. This will assist me in making sure all the new and modified streets are updated in OpenStreetMap.
The street segments should have a unique identifier, so my only thought is to write a Python script which opens both Shapefiles, finds any segments in the new that aren't present in the old (added segments), and the reverse (removed segments), and then loop over matched pairs comparing their constituent coordinates to see if anything has changed.
Is there a way to do this using QGIS?


texts: ArcGIS Random Point generator producing same points every iteration

texts:
I am generating random points to sample within a stratum.
Data Management Tools > Feature Class > Generate Random Points
This works just fine, except that I noticed the points I generated for sampling this year were exactly the same as last year. Every time I generate a new set of points, they are the same points each time.

Is there a setting that will fix this?

I am using ArcGIS Desktop 10.3.1.


texts: 3D interpolation of pollutants across water depth

texts:
I have point data representing tissue sampling performed by NOAA NCCOS for Gulf of Maine waters. I am assuming a temporal resolution of one sampling point per location for one year. I have searched for literature to address how to interpolate this data across water depth, but have yet to come to a conclusion.

A simple kriging shows patterns of overlay, which as a simplified approach, is fine. Are there more refined methods to determine the distribution of pollutants using marine tissue samples? I would prefer methods that have been published or applied.

I am using ArcGIS 10.2 and I am open to any open source software. I am also proficient in Python and ArcObjects.


texts: Error when using input() pyqgis (RuntimeError: input(): lost sys.stdin) QGIS 3

texts:
I need the user to input a value. I'm using input() on qgis python console, but getting this error:

  RuntimeError: input(): lost sys.stdin

How can I solve it?


texts: Exporting jpeg around polygon boundary using Python?

texts:
I have a python code that goes through my map document and exports a jpg for each polygon grid.  Right now it clips the image to the boundary of the each polygon.

Does anyone know how I can make it clip a small amount of data around the grid?

I only want a 1-2% more added but I'm not sure if this can be added to my existing code or do I have to create new polygons that are all a small amount bigger so I can capture this extra data.

Below is the code I am using right now.

import arcpy
... mxd = arcpy.mapping.MapDocument(r"C:\Users\markd\OneDrive\Documents\GIS\PalmSprings\PalmSprings_TileExport.mxd")
... df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
... lyr = arcpy.mapping.ListLayers(mxd,"L8_PS")[0]  ##shapefile layer name called grid
... rows = arcpy.SearchCursor(lyr)
... x=1
... for row in rows:
...     df.extent = row.Shape.extent
...     arcpy.mapping.ExportToJPEG(mxd, r"C:\GIS\PalmSprings\L8_PS\ " + str(row.getValue("PageName")) + ".jpg", df, df_export_width = 1200, df_export_height = 1200, resolution = 96, world_file=True, jpeg_quality = 100)
...     print x
...     x+=1
... del mxd


texts: Programmatically associating layer with custom UI using PyQGIS?

texts:
I know it is possible to associate a layer with a custom UI for editing feature attributes using the QGIS interface. But is this possible with Python?

This link shows how to achieve the result manually. But I would like to be able to create layers automatically and have them use custom UI without having to redefine the ui file and method each time.


texts: How to digitize a "ringed" topographic-style choropleth map?

texts:
I have a choropleth map of plankton distribution in the ocean that I am digitizing. However, as you can see, there are many features within features creating rings, almost like mountains on a topographic map. I'm using QGIS to digitize this map.

At the moment my game plan for digitizing these is to:

Digitize the innermost rings using Freehand editing plugin - "> 500" according to the scale
Digitize everything surrounding those - "250-500" according to the scale
Clip all "> 500" patches out from the "250-500" shapefile

So on and so forth. This does mean I'd have many shapefile layers and much clipping to do, though it seems faster and less error-prone than simply tracing out each feature and contour line by hand.

My question is - before I commit too much time and effort to this clipping technique, does anybody in the community know of a better or faster way of doing this? Or is simply manually tracing each contour line the best way to go?

I'm worried that I may be reinventing the wheel,  trying to figure out a problem that may be common knowledge to somebody who has much digitizing experience.

PS - if there is a better term for this map than "psuedo topographical choropleth", and I assume there is, please update my post/let me know! I don't know what this style of map is called.


texts:
I want to find villages/cities which are ends. That means: I need this villages which are linked to maximum one village on a distance of X km.

I'm just wondering how is this algorithm is called (I guess someone created something similar before).

Is there a tool for this? How is this concept called?

E.g.

 City ------------ Village 1 ------- Village 2 ------ Village 3 - - - [ Mountain]
                    \                 \                                \
                     Village 4         |                                Road
                     |                  \                            in mountains
                     Village 5----------Village6---------Village7 - - - - -

As long the Road in mountains is longer than Xkm, Village 3 and Village 7 would be considered by me terminal, because they are linked to one village (V2, respectively V6).

Otherwise, I will want to build one with the open data we have.

What I've tried:

Downloading the OSM data for my country (Romania) and importing the villages and cities into the database. Using the geolocation functions from the database I'm able to find the villages which do not have more than X villages in the radius of R km.

However this is not a solution for my case because in my cases a village may be on the other side of the mountain, like in the example above, but there is no good way to it (or no way at all).


texts: Geocode service findAddressCandidates returns no candidates with 10.0 address locator

texts:
I posted this on an esri forum but nobody touched it. We've recently updated a geocoding service and I'm now unable to return addresses to my Java application. The old service is built with a 9.3 address locator and is served from a 10.0 ArcGIS Server.

The following code is used to return intersections after the user enters street names.

String url =
MessageFormat.format("{0}/findAddressCandidates?Street={1}%7C{2}&City=&Single+Line+Input=&outFields=&outSR=&f=pjson", geocodeServiceUrl, street, crossStreet);
GetMethod httpGet = new GetMethod(url);
int statusCode = client.executeMethod(httpGet);

The following JSON is returned with the old service:

response={
"spatialReference" : {
"wkid" : 2278
},
"candidates" : [
{
"address" : "EAGER ST | CALLAGHAN AVE",
"location" : {
"x" : 2132338.7499425262,
"y" : 13698061.000045776
},
"score" : 69,
"attributes" : {

}
}
]
}|

However, when I switch to the new service no candidates are returned in the response:

[#|2013-01-09T13:34:51.176-0600|INFO|sun-appserver2.1|null|_ThreadID=18;_ThreadName=httpSSLWorkerThread-8080-3;|geocodeIntersectionUsingREST: url=same URL/ArcGIS/rest/services/Locators/PWStreets/GeocodeServer/findAddressCandidates?Street=eager%7Ccallaghan&City=&Single+Line+Input=&outFields=&outSR=&f=pjson|#]

[#|2013-01-09T13:34:51.959-0600|INFO|sun-appserver2.1|null|_ThreadID=18;_ThreadName=httpSSLWorkerThread-8080-3;|geocodeIntersectionUsingREST: response={
"spatialReference" : {
"wkid" : 2278
},
"candidates" : [

]
}|#]

Strangely, when we run the same query using the ArcGIS Services directory the correct addresses are returned in JSON. The requests from the app and the services directory are identical. So it seems that only in my Java application and using a 10.0 service and address locator are no results returned.

Until I figure this out we are stuck using the old service built with a 9.3 address locator. Very confused.

Thanks!