texts:
I work for a local government and am in the process of cleaning up / updating addressing data. I have a MXD for the purpose of 'point addressing' and 'centerline addressing'.

It is the centerline data that I am concerned about. My workflow can either consists of "rolling my sleeves up" and manually digitizing the data or what I wan to do is find an automated way to complete this work.

The city has theoretical and actual address ranges. The issue is that the actual address range data is filled in the theoretical fields, thus the theoretical address ranges are wrong.

What I am trying to accomplish is writing or using attribute assistant to identify all actual address records that are Null and copy the records from the theoretical address fields to the actual address fields. Further since the theoretical address fields are partially wrong I would like to have a TRIGGER event that if the process where the Null values are converted = TRUE then the last two digit of the Left_to and Right_to fields for the theoretical addresses will be changed to -98 and 99.

Example:

Actual Address Range:
Left side of road = 800-832 ; Right side of road = 801-835

Theoretical Address Range:

Left side of road = 800-898 ; Right side of road = 801-899


texts: How to calculate the sum of all cells in a raster layer upstream of a specific point?

texts: I have a rainfall map which covers large area. However, I'm interested in the total value of rainfall over a specific subcatchment. If I don't have the shapefile which defines the boundary of this subcatchment, is it possible to calculate the sum of the rainfall values upstream of a specific point in ArcGIS? or Do I have to delineate the catchment upstream of this point first?


texts: Dissolve using ogr2ogr gives 'no such function: ST_union'

texts:
I am trying to dissolve a shapefile using ogr2ogr but I am getting error no such function: ST_union:

Command:

  C:\ms4w>ogr2ogr "C:\output\test-dissolved.sh p" "C:\input\test.shp"
  -dialect sqlite -sql "SELECT ST_union(geometry) FROM test "

Error:

  ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_union(geometry) FROM test ):   no such function: ST_union


texts:
I have shapefile which is too large to be served as vector. Also, I need it only as background tile, not as layer with attributes.

I thought of creating tiles in MapTiler as I already bought it. I see MapTiles accepts raster input. Where should I create rasters from shapefile? I need it for only 3-4 zoom levels. What would be next step? I have a very good server, should I uploaded it to my server? Should I process it via Geoserver as WMS? What would be best approach if I need it in Leaflet application?


texts: QGIS 3x capturing coordinates for multitude of objects at once

texts:
I would like to ask about the possibility of capturing coordinates for more than 1 object in QGIS.
I used to use the coordinate capture plugin, but what is the problem there I have to collect coordinates one by one.
Is there any way to collect coordinates i.e. for all points from the given layer and export them to excel? Following their names or different features from data attribute table?

I will by more than happy if someone clarifies it for me.
UPDATE:
Following the answer below, I tried to run the $x expression, but it returns NULL.

What am I missing here?


texts:
How do I set the visibility scale for a line KMZ created from ArcMap (Layer-to-KMZ)?

I'm sure I would have to go into the KML code and add, but add what.  That is to say as I zoom in and out using Google Earth I want certain layers to be visible and other not based upon appropriate level of visible details.


texts: Conflicts while installing postgresql9.3-postgis2.1,QGIS and gdal-bin

texts:
I want to use postgresql9.3-postgis2.1 and gdal-bin parallel .

I have refer instruction here for installing postgresql-9.3 on ubuntu 12.04.
I have install qgis from here

While installing gdal-bin it will remove postgis. As I explore I found similar question which says ubuntugis-unstable repository and postgresql repository will not work together.

Then how can install postgresql9.3-postgis2.1 and gdal-bin parallel. Can someone help to solve the problem , Any suggestion will be great help for me.

Thanks


texts: Zonal Statistics as table in model with two iterators using ArcGIS Desktop?

texts:
I am using ArcMap 10.4 (with Spatial Analyst.)

I have a feature layer with ~11.000 partly overlapping but similar features. For each feature I want to extract some values from ~25 underlying rasters.

I set up a model that iterates the raster folder and applies the needed values to the features. It works as intended, except for Zonal Statistics as a table doesn't handle overlapping features.

I have tried the following in order of solving this problem:

Set up the iterators in submodels, but they only run once. In order of fixing that I tried making them preconditions for the forthcoming tools. That didn't solve produce intended result.
I tried simply deleting every second feature, by MOD(layerFID,2)=1, which was close to correct, but isn't really a solution.
I set up an iterator that handles the overlapping features, by splitting the into separate shapefiles, performs Zonal Statistics as a table on them, and merges them again. It works, despite of being very slow. And then I can't iterate the raster folder.
I tried running Zonal Statistics 2 (and 3), with and without "arcpy.env.overwriteOutput = False". I yields odd results. As far as I know COUNT should be the same for all the zones when the area and the raster is the same.

I don't see any decent solutions that doesn't involve two iterators.


texts:
Even though I added the Arc Hydro Tools to my toolboxes I cannot open and access the Arc Hydro Tools.

[


texts: How to enable buttons of a value table in python toolbox?

texts: Recently I asked a question in stack (How to search on values of a value list in python toolbox?). This question is about enable some default buttons in the tool. I wonder why the add, delete and order buttons not working for me. I want to order entries(rows) by order buttons but the order buttons only select records. Is that any config to enable these buttons? The full image of this tool can be found in the past question.


texts:
In the advanced digitizing toolbar, QGIS has these 3 tools: Split Features, Split Parts, and Merge Features. Is there a tool or equivalent algorithm for the logical fourth one: "Merge Parts"?

This would have the opposite effect of Split Parts, which takes one part of a multipart feature and splits it into two adjoining parts of the same feature. Conversely, "Merge Parts" would take two (or more) adjoining parts of the same feature and make them a single part of that feature. It would not affect the singlepart/multipart status of the feature (or the layer as a whole).

In the image below there are two multipolygon features, which I realize I should have called Feature A and Feature B (rather than Polygon A and Polygon B). Feature A has 2 parts, and Feature B has 4 parts; note that parts 1 and 2 in Feature B share a common boundary. I've selected Feature B, and I'd like to execute a tool to join parts 1 and 2 into a single part, as shown in the "after" section. Feature B would then still be a multipolygon, but would have only 3 parts instead of 4.

The closest I can find is the "merge lines" algorithm; however this would not work with multipolygon features. It also merges all connected parts of a line feature, instead of just selected ones. (Is it even possible to select individual parts?) Is there anything else I'm missing?


texts: Removing Layer on ArcGIS API For JavaScript>

texts:
I am adding a markers layer called layer1 like this to map

function drawPoints(mapInfo) {
  layer1 = new esri.layers.GraphicsLayer();
    for (var i = 0; i < mapInfo.length; i++) {
           var projects = mapInfo[i];
           var project = new esri.geometry.Point(projects.Longitude, projects.Latitude);
           project = esri.geometry.geographicToWebMercator(project);
           var symbol = new esri.symbol.PictureMarkerSymbol("img/map/marker.png", 18, 18);
           projectInfoTemplate = new InfoTemplate();
           projectInfoTemplate.setTitle("Project Details");
           projectInfoTemplate.setContent('

 ');
           var projectsG = new esri.Graphic(project, symbol).setInfoTemplate(projectInfoTemplate);
           layer1.add(projectsG);
           }
          map.addLayer(layer1);
 }

now in next request I need to clear map so I used the

map.removeLayer(layer1);

but this is causing error because the layer1 still not created at first request. Now I need to check IF the map has a layer called layer1 then removeit. Here is a pseudo code of what I need to do:

if(map.has/contains/include(layer1){
  map.removeLayer(layer1);
}

can you please let me know how to do that?


texts: MultiRing buffer plugin to Python script in QGIS

texts:
Is there a way to get the MultiRing buffer plugin into a Python script?
I want to use it in the QGIS Modeler, but I cannot use plugins in there, or to get multiple buffers with varied (not equal) steps into the Modeler.


texts: Converting projected coordinate system worldfile to geographic coordinate system worldfile?

texts:
I have an image with the worldfile like this:

0.298582141739
0.000000000000
0.000000000000
-0.298582141739
1283836.327077804830
6134835.890168172310

I think it is the projected coordinate system. I want to convert it to geographic coordinate system, because I want to know the latitude/longitude of each pixel of the image.
How to convert that?
I check the worldfileread in the matlab, it can read both. But I did not find how to convert.

Could anyone know about that?


texts: How to display together a shapefile (country basemap) and points (from CSV) via geopandas?

texts:
So, I wrote this code that tries to do it but I end up with the map and points being displayed separately

import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame
from shapely.geometry import Point
import fiona
import matplotlib.pyplot as plt

def loadPoint():
    df = pd.read_csv(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\emscPhilippines2008to2020_lonLat.csv')
    geometry = [Point(xy) for xy in zip(df.Longitude, df.Latitude)]
    crs = {'init': 'epsg:32651'}
    geo_df = GeoDataFrame(df, crs=crs, geometry=geometry)

    geo_df.to_file(driver='ESRI Shapefile', filename='eqpts.shp')

def plotPoint():
    map = gpd.read_file(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\PH_provs.shp')
    map.plot(cmap = 'jet', column = 'NAME_1', figsize=(10,10))

    pt = gpd.read_file(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\eqpts.shp')
    pt.plot()

loadPoint()
plotPoint()


texts:
How do I do this?

Columns -    Sitecode       Acres
             271730-001     .75
             271730-002     .50

I'm trying to label objects by multiplying the Acres column by "2.5".
As an example using first row from above:

????(Acres*2.5, ??)  should label as "1.875"


texts: How to join/relate points to polygons by field

texts:
I am using ArcMap with a basic license. I have a layer of polygon parcels and a layer of point turnouts. The turnouts deliver water to the parcel, but are not located within the parcel.
My parcels have a text field for the Turnout Unique ID, for example (2019,2020,3410). I spoke to an ESRI rep and they mentioned a relationship class tool that is available under a standard license. I want to be able to click on a parcel and the associated turnouts are shown. Is that the only tool to do this task?


texts:
I've got my final raster tile map (via vrt) into tilemill and am trying different mbtiles export options. No matter what boundary I specify for export - I always end up with the beige top and bottom bands.

How do I get rid of them?


texts: New field added whenever I export shapefile to the geodatabase

texts:
I don't know why this is keep happening, but I created shapefile in QGIS and opened the shapefile from ArcGIS. I needed to create geodatabase to put in all of my shapefile, but whenever I do that, point shapefiles are fine, but my line shapefiles automatically create new field.

The main problem is that I cannot delete that new field in my attribute table, but I want them to be disappeared.

It is just one field that is keep showing up. I got the shape_leng for the shapefile's length(m) but when I export that into geodatabase, the new field shape_length comes up. It doesn't let me delete this field and values are very different.

Is there any ways to fix this or delete this field permanently?


texts: How do I remove a Reader from my Google Earth Engine?

texts: I accepted the Repo shown in the Reader category.  How do I remove it?


texts: Getting zonal statistics of polygons not touching center raster cells in PostGIS?

texts:
I'm working a postgis processing zonal stats from raster and I couldn't manage the problem of raster with so low resolution (like 9km width) and small polygons.

My problem is, these rasters are big (for all South America) if I try to resample the pixels to smaller ones, it will spend a lot of time.

There are a possible solution using the ST_SummaryStats() style, like my example code bellow.

I try this code and it only get stats when the polygon covers the centroid of pixel image.

SELECT gid, St_SummaryStats(St_Union(ST_Clip(ST_SetSRID(rast, 4326), 1, ST_SetSRID(geom, 4326), 0.0, true))) FROM raster_image, shape_polygons WHERE st_intersects(ST_SetSRID(rast, 4326), ST_SetSRID(the_geom, 4326)) GROUP BY gid;

The Example:

The results of Stats


texts: How to load spatial data into ArcGIS Javascript API application from SQL Server 2008?

texts: I have geometry data in SQL Server 2008 R2. How do I load that spatial data using the ArcGIS Javascript API?


texts: Georeferencing PNG with transparency via gdal produce "green pixels" in QGIS?

texts:
This is just intriguing while not blocking my way, but I wonder what is going on and would like to get feedbacks from more experienced users.

1) Given a trans.png file (1.1MB) looks good and conserve transparencies (Stackexchanges displays below an altered version as .jpg) :

2) I restore georeferencing using :

gdal_translate -a_ullr 67.0 37.5 99.0 27.0 -co COMPRESS=LZW  -co ALPHA=YES       ./trans.png ./trans.gis.tif

3) The output trans.gis.tif (1.6MB) looks good and conserve transparencies (Stackexchanges displays below an altered version as .jpg):

4) Loading it in Quantum GIS, it get well positionned (georeferencing works!), but the image is ill colored with the darkest areas (and only them!) becoming green :

Is this green artifact normal ? How to fix it (so I may see the expected black pixels) ?

Edit: file added in 1)


texts: Is spatial reference information suceptible to be translated?

texts:
In an English application that uses ArcObjects API from ESRI, the information regards the spatial reference system is displayed in English as the API gives it.

If the application is translated to another language (let's say Russian), the information regards the spatial reference should be translated too? I mean, terms like "Mercator", "Degrees", "Magna", "WGS1984UTM30N" should be translated?

Does any one know if this information appear in Cyrillic alphabet in ArcGIS in Russian? I wonder, Is there a Russian version of ArcGIS Desktop? if the answer is positive, how does it show the information of spatial reference systems?


texts:
I am looking for a source for data to determine voting precinct boundaries in the United States. My thought is that they might be on the Census website but under a name that doesn't make it obvious.

Does anyone have a working link (or set of instructions) that will allow me to find that information?


texts: Running Random Forest Classification on WorldViev-2 image in R - problem

texts:
I have a problem with performing RF in R for a WorldView-2 Scene. I've never ever before used R for remote sensing classification, so I simply followed what is written here.
I've loaded scene (8-band TIFF) with raster::brick and called it abcd, loaded simple 4-point, 2-class shapefile with OGR, run randomForest model and called it wv2rf and tried to use predict:

predict(abcd, wv2rf, "rftes.img", index=1, na.rm=TRUE)

R returned an error, that I cannot handle with:

Error in `[.data.frame`(blockvals, , f[j]) : undefined columns selected

I guessed that it has to do something with index value, but I tried to change it and nothing worked

whole code:

#loading data
abcd <- brick("13AUG20094646-M2AS-13EUSI-1283-01.tif")
sdata <- readOGR(dsn="mypath\\poligony testowe", layer="training")
v <- extract(x=abcd, y=sdata, df=TRUE)
sdata@data = data.frame(sdata@data, v[match(rownames(sdata@data), rownames(v)),]
sdata@data[3] <- NULL
#classification
wv2rf <- randomForest(x=sdata@data, y=(sdata@data[,"klasa"]), ntree=10, importance=TRUE)
predict(abcd, wv2rf, "rftes.img", index=1, na.rm=TRUE)

Here are files with my sdata, wv2rf and abcd variables, should work with load("filename") in R. Hope they help you help me.

http://speedy.sh/9f7cf/rf.zip


texts: How to loop in PostGIS to add the next lat long to a route

texts:
I want to loop over the results of a query and then pass row(n)'s longitude, latitude and row(n+1)'s longitude and latitude to another function.

I have a query which returns a set of coordinates from a way-point table with no geometry field.

SELECT longitude, latitude FROM way-points ORDER BY date

I am also using the pgRouting extension and I have replicated the pgr_fromAtoB function to allow me to generate a route between two coordinates.

Essentially, I have a table of waypoints (GPS lat/lon) and a separate table of roads. I need to get each GPS lat/lon, then get the road (or roads) that make up the route to the next waypoint.

My idea is for each waypoint pair (row(n) origin, row(n+1) destination) in the waypoint query, call the pg_Routing function to find the nearest road vertices's (NNS), then route between the two.

All I have to do is pass the input: x1,y1,x2,y2 and I get back the list of geometries.

SELECT geom FROM pgr_fromAtoB('road_table',-122.662,45.528,-122.684,45.514);

Ideally I would like to have this automated so that instead of specifying the values of x1,y1,x2,y2 I can just use a SELECT query that returns the coordinates.

So maybe something like:

BEGIN
-- Generate x1, y1, x2, y2

-- Set origin and remember start point
    x1 := 0;
    y1 := 0;

    FOR rec IN EXECUTE 'SELECT * FROM waypoints WHERE user_id = 5 ORDER BY waypoint_timestamp'
    LOOP
    -- first iteration sets x1, y1
    IF ( x1 = 0 AND y1 = 0 ) THEN
        x1 := rec.longitude;
        y1 := rec.latitude;
        geom1 := rec.geom;
    ELSE
        -- If the origin and destination geom are the same don't route
        IF geom1 = rec.geom THEN

        ELSE
            -- geoms are different so furuther iterations set x2, y2
            x2 := rec.longitude;
            y2 := rec.latitude;
            -- call the pgr_fromAtoB function (road table hard coded)
            EXECUTE 'SELECT ST_Union(f.geom) FROM pgr_fromatob(''roads''' || ',' || x1 || ',' || y1 || ',' || x2 || ',' || y2 || ') AS f' INTO geom;
            -- now set our new start point equal to the previous destination
            x1 := x2;
            y1 := y2;
        END IF;

        RETURN NEXT;
    END IF;
    END LOOP;
    RETURN;
END;

Do you know how to write loops in SQL and can help me refactor this into something that actually performs well in PostGIS?

Perhaps you know a faster way to do this without loops. Performance is my biggest consideration.


texts: QGIS 3D and QGIS Webmap, Javascript; Search Bar

texts:
I exported my webmap, using QGIS2web.

The QGIS2web already has a search bar, but the placement icon when I search, as it applies to my project, always shows below my uploaded raster images. When I zoom in it vanishes, because the size of the search icon reduces, and lies below my raster files. I used a series of images to explain.

Figure 1: I highlighted my uploaded polygon at the right bottom corner. The underlying map is OSM.

Figure 2: I searched for a feature, and the box suggested related features for me in the top left area.

Figure 3: haven clicked the enter icon to validate my search, I unchecked the images, showing the underlying OSM. You should be able to see the blue dot representing the response of the search

Figure 4: As I zoomed inward, the blue dot becomes clearer, radiating from beneath my polygon that I didn't uncheck.

Figure 5: As I zoomed further inward, the search response icon becomes smaller than my individual raster, therefore unable to show, as its default is to show under my polygon (this is what I want to correct. I want it to show above my images.)

Figure 6: To prove that it was underneath, as I unchecked my raster images, it became more clearly visible.
Isn't it possible to just have the search function highlight my search feature in a different color, instead of using a dot to point to it, as is the case here?


texts: Function in R to get difference between two sets of polygons comparable in speed to QGIS

texts:
I have two very large polygon shapefiles for which I need to calculate the difference between the two shapefiles. One is a shapefile of forest polygons, the other is a shapefile of road polygons. I want to remove the area of the forest polygons that intersects with road polygons, while maintaining each polygon as a separate feature and its data.

After much trial and error, I found that erase using the sp format or using st_erase using the sf format performs the operation that I need.

st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))

However, the processing time is really long. With a small sample of my dataset, I have processing times several minutes long. (And R crashes with a tenth of the dataset). Meanwhile, the function Vector-->Geoprocessing tools--> Difference in QGIS takes 83 seconds with a tenth of the dataset.

erase(forest_sp, roads_sp)
#2.83488 minutes
forest_sp-roads_sp
#2.8698 minutes
st_erase(forest, roads)
#58.77sec

Is there a function in R performs this operation quickly? Or is R inferior to QGIS in processing times for this kind of spatial function?

Background

There are multiple questions on gis.stackexchange that deal with this type of problem. Some are with the sp format and suggest gDifference or erase (or its analog polygon1-polygon2):

Erase, gDifference

Clip polygon and retain data?

gDifference

Remove parts of line that fall within a polygon

Identify area within one shapefile not in another shapefile in R

Reverse clipping (erasing ) in R?

Other questions suggest solutions with the sf format:
st_difference

Reverse clipping (erasing ) in R?

Using simple features (sf) in R, how do I erase polygons overlapping with another layer

ogr2ogr equivalent of QGIS Union

All state that they perform that same thing (remove intersecting portions) and the official tag definition on gis.stackexchange states that a difference operation removes sections. However, I have found that often they don't.
For example, gDifference with byID=FALSE that merges all into one MULTIPOLYGON does remove overlapping sections but gDifference with byID=TRUE that keeps polygons separate includes overlapping sections. With sf, st_difference also keeping overlapping sections and ends up with 10x as many polygons as my original file.

However, erase and the helper function st_erase do work as expected.

Difference

  A "difference" operation is a geoprocessing operation that completely
  removes one overlapping polygon and the intersection of that polygon
  from the other. It uses one polygon to "take a bite" out of another.

Erase

  The action of removing data attributes or features through the
  processes of positive or negative selection and deletion.


texts: Selecting set of vertices (not complete feature) and pasting into another feature in same or different layer using QGIS?

texts: I want to select a few vertices in QGIS 3.2.2 (for example a bend in a river) and copy that into another feature in the same or in a different layer.  The pasted vertices would be added to existing vertices in a feature.


texts:
How can I add or subtract from a date in Google Earth Engine?
See my codes here-

var imagedate = ee.Date(collectionLT05_1991_1.get('system:time_start'));
var imaagedate1 = imagedate+1;
print(imagedate1);

I want the day after the "imagedate".
Is there any specific library for this ? or suggest me the appropriate codes.


texts:
The HTML annotations in qgis2 seem to fit the bill of what I'd like to do: have a full-custom data-based drawing in a callout.  Problem: I can't find documentation anywhere on how to write the html file to incorporate placeholders to take data values from the annotated object.  Web searching shows a few discussions that say something about placeholders, which is the only evidence so far that it could actually do the trick.  Can someone please point me to documentation on writing the html file used by html annotations?  Maybe I'm just missing something obvious?

Thanks


texts: Filter NDVI polygon region Google Earth Engine

texts:
I'm trying to get the NDVI specs from a polygon region around SF and the Bay Area but the image isn't working. It is being displayed in a region next to the polygon or just comes out black
My code is as follows and the variable "geometry' is the polygon I have drawn on the map

// Get the least cloudy image in 2016-17.
var image = ee.Image(
  l8.filterBounds(geometry)
    .filterDate('2016-01-01', '2017-12-31')
    .sort('CLOUD_COVER')
    .first()
);

// Compute the Normalized Difference Vegetation Index (NDVI).
var nir = image.select('B5');
var red = image.select('B4');
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');

// Display the result.
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
Map.addLayer(ndvi, ndviParams, 'NDVI image');

Export.image.toDrive({
  image: ndvi,
  description: 'NDVI_CA',
  scale: 30,
  region: geometry
});```


texts: Calculate and add area of overlapping features within polygons

texts:
I have a polygon layer of watershed boundaries

And a single multipart feature that contains all public lands for the state.

I need to find the area of public lands within each watershed and add it as a field to each of the watershed polygons.

For example, for the Black river, I just want the area of public lands (tan) that is contained within the polygon.

I assumed I would use a spatial join from this, but I can't figure out how to get the area within each watershed as opposed to the area of the giant multipart public land feature.


texts: OpenLayers, Geoserver, Postgis, GITHUB and content management

texts:
Im not sure if this is the right forum for this question.

I am working on a project to look into developing a prototype mapping applications for the agriculture industry. I am using asp classic, openlayers, geoserver and postgres for the application atm I plan on changing from asp classic to php soon.

Could  a CMS be used to manage page layouts, map layers styles, user profiles, user permissions etc in my current aplication.

Can GIT repository be used to store code, and manage code versioning used in the application?

I have no idea where to start. Has anyone got any advise or links to share. I have done some searching on the web but there is just some much to trawl through.
I do have both cartaro and drupal installed, but I can find any simple way to use them for my application. Here is a linkhttp://54.252.203.210/FARMWEB/FARMWEB2.0_DEMO_STEP3_PG.asp?farmid=354 to what I have developed so far.


texts: QGIS - point to raster - Summing attributes as value of raster

texts:
I'm trying to produce a raster layer. The aim is to show the summed heat consumption of every cell (100*100m).

I created an point layer with all buildings and an attribute with their heat consumption. Using the Rasterize (Polygon to Raster) funtion I can choose the attribute field heat consumption, but the output doesn't show the sum of the attribute heat consumption in every cell. The value of every cell is a random value of one building inside and not the sum of all buildings inside.The problem is I need the sum of every building inside the cell. I tested a lot of different ways to find a solution, since I'm still on beginner level, it's hard for me to find out a way to calculate this.


texts: Adding raster layer with RPC file in ArcGIS for Desktop for ortho rectification with DEM file?

texts:
I have a Pleiades image with its RPC file with .xml extension and a DEM file of the same area, but when I import them the orthorectification button is not active in display tab of property window.

How can I add them?


texts: How to extract two values from an 8bit raster and export to new raster with FME 2015?

texts:
I have an 8 bit raster from which I want to create a new raster which just contains the cells containing the values 20 and 24 of the original raster.

I have SAFE FME 2015 but am unsure which transformer to use.


texts: Search for untagged ways in overpass-turbo

texts: Is there  a way to search for untagged ways that are not part of a multi-polygon relation using overpass-turbo.eu?  I am trying to locate all ways in my county that were not tagged (i.e. untagged highways, untagged buildings, pools, etc.)


texts: pyqgis splitGeometry practical working example

texts: I want to get splitGeometry working example in qgis where i can find it. Official documentation have define here but this is not enough. Where i can get its working example in order to understand it and use in my project.


texts: Calculating shortest path between each lat/lon point in one table and each lan/lon point in another in PostgreSQL?

texts:
I have two tables, both of them have only 2 columns. The first one contain only latitude and longitude coordinates of my source points and the other one only latitude and longitude coordinates of my target points. For each source point I have to find the closest (along the streets) target point (and tell what's the distance). I have already developed a function give_me_cost, that finds the closest node and performs the routing:

create or replace function nearest_node (double precision, double precision) returns

integer as
$$
SELECT id::integer as source_id FROM roads_vertices_pgr
                        ORDER BY the_geom <-> ST_GeometryFromText('POINT('
                        ||$1|| ' ' ||$2||')',4326) LIMIT 1;
$$ language sql;

CREATE OR REPLACE FUNCTION give_me_cost(double precision, double precision, double precision, double precision)
  RETURNS decimal(8,8) AS
$BODY$
with subset as (
SELECT seq, id1 AS node, id2 AS edge, cost
  FROM pgr_dijkstra(
    'SELECT gid as id, source, target, st_length(geom) as cost FROM roads',
    nearest_node($1,$2),nearest_node($3,$4), FALSE, False
  ) as di
)
select sum(cost)::decimal (8,8) from subset

$BODY$
  LANGUAGE sql VOLATILE;

Now I have to use it in some kind of loop, but don't know how. The most preferable solution would be the function, where I can input my source table and my target table and get an output of a table with sources, the closest targets and the distances between them.


texts:
I'm writing a C++ application using QGIS API and I would like to improve the performance.

I want to render 5 layers:

one raster layer using gdal to load Bing/Google Map tiles
4 custom QgsMapLayer (using plugin layer type) to draw different lines/polygons/points in WGS84.

For now, each of these layer is added to the MapCanvas, using:

QgsMapLayerRegistry::instance()->addMapLayer(customLayer);
QList <QgsMapCanvasLayer> myLayerSet;
myLayerSet.append(QgsMapCanvasLayer(cutomLayer,TRUE));
setLayerSet(myLayerSet);

On the canvas, I can have map tools, like panning, zoom in/out.
I use custom layers because I want to easily decorate lines and polygons, using QPainter.
This works, and I'm able to display what i want, where i want (using coordinate transforms from WGS84 to pseudo mercator).

My problem is: one of my custom layers draws a GPS point with some circles around it.
This layer is refreshed every second.  To refresh my layer, I use the canvas->refresh(), but that refreshes all other layers as well, and by consequence, the GDAL layer. This takes a lot of time and it's not necessary.

I have read some posts here mentioning the use of layer->setImageCache(), or by sub-classing QgsMapCanvasItem. Unfortunately, none of theses solutions works as I expect. setImageCache apparently doesn't work at all, or I don't know how to use it properly.

Sub-classing QgsMapCanvasItem allows me to implement the layers that need frequent refresh, and refresh them without refresh Raster Layer but:

when panning the map, drawn objects don't follow the map (as a rubberBand can).
when i use zoom in/out, the selection rectangle erase what was drawn on my canvasItem, without triggering a redraw.

I don't see how to do that, even if this subject seems relatively common.

Am i on the right way or is there a better way?

EDIT :

thanks to Nathan, my problem is resolved, and i confirm that custom QgsMapCanvasItem is the solution to independently refresh some layers.
just keep in mind to use toCanvasCoordinates(), and implement boundingRect() properly.


texts: ArcMap bookmarks for legends, scale bar etc

texts:
I am exporting several GIS dataset maps from a set of several files I have. My images show always the same geographic area. My problem is that, moving from one image to the other one, the legend, north arrow, scale bar and so on are always located in different location (e.g., one is slightly more close to the upper right bottom than the other one).

Is there a way to define a "bookmark" in ArcMap to place these elements always in the same position of the paper image in layout view?


texts: Python Add-in- Message Box to display area of selected polygons

texts: I am looking to create a Python AddIn which calculates the area of selected polygons.  I would like to be able to select the features, click a button on the toolbar, and a Message box displays the sum area.  Is this possible without writing the area to the table?  I hope that makes sense.  I am very new to Python.


texts: Is it possible to show multiple infowindow?

texts:
I just started working with API Arcgis.
I can not find a way to leave open the popup that I add with:

mymap.infoWindow.setContent();

I wish that every click of a POI then opened his infoWindow.
Is there a way? Documentation?


texts: ArcPy Data Access Insert Cursor not writing geometries to feature class?

texts:
I'm using PyScripter and running ArcGIS 10.2.1.

I am having troubles with the Insert Cursor in the arcpy module.

I have an empty feature class in my geodatabase and have created the short script below to try and find the issue, but have not been able to. When I run the script, it executes fine, but when I look in the attribute table of the feature class afterwards it is still empty.

import arcpy

fc = r"C:\Users\djh\Desktop\topo_map\test.gdb\well_location"
cursor = arcpy.da.InsertCursor(fc, ["SHAPE@XY"])
xy = (206901.75, 5997594.47)

cursor.insertRow([xy])

I have copied this script straight from the "writing geometries" page on the ArcGIS help website, so I don't think it's a problem with syntax. Please let me know if there's a way to get around this issue.


texts: Installing Python modules for QGIS 3 on Mac

texts:
I want to install some third-party Python modules for QGIS 3 on Mac, however, I have no clue how to do it, since I am not so much familiar with different Python environments.
Is there any similar thing to the OSGeo4W Shell on Mac, or is it possible to do via the terminal?


texts:
Is there a crosswalk table from USDA agricultural district to FIPS code?

USDA agricultural district https://www.nass.usda.gov/Charts_and_Maps/Crops_County/boundary_maps/asdsbig.gif

I have checked https://www.agcensus.usda.gov/Publications/2012/Online_Resources/Ag_Atlas_Maps/mapfiles/ag_co_metadata_faq_12.html#what

It seems that I have to do some complex calculation to get the crosswalk table. I do believe there is a crosswalk table from USDA agricultural district to FIPS code.


texts: How do I make a variogram of my model residuals in R?

texts:
I have a csv file containing LST values for several study sites, canopy cover for each study site, and impervious cover for each study site, and the geographic x/y coordinates for each study site. I have created a linear model as: LST ~ canopy cover + impervious cover.
I want to check for spatial autocorrelation by plotting a variogram of the model residuals. I've been trying to use the variogram() function in gstat but can't quite figure out the proper way to write the code (in particular for the locations argument).


texts: Store global variables in ArcPad that store user information

texts:
I know that it's possible to store global variables in ArcPad from reading the help files.  I have several field staff that are out everyday -- generally we have them select their name from a "Specialist" domain for each data point they create.

Is it possible to store a table somewhere in the background of ArcPad that stores the user's information that can be used to auto-populate fields such as Specialist (user-name)?


texts: Splitting polygons layer by using line layer not working in QGIS?

texts:
I am using qgis in project and trying to split polygons layer by using a line layer.
I used the the splite tool in qgis as explained in many posts.
The output after splitting the polygon layer is deleting all the left side polygons from the  splitter and kept the polygons at the right side. One side kept and the other side disappeared.

How to keep both polygons sides from the line splitter?


texts:
We just installed portal 10.6.1, completely new environment... almost everything appears to be functioning normally except point symbology configuration will not load within a webmap

.

We're seeing the same behavior in the default symbology settings within the item settings for point layers.  I also checked out the REST endpoint and noticed that the drawing info section looks different that in previous versions... but perhaps this is just a new syntax?

.

Polygon and polyline symbology appears to be working fine--seems isolated to point layers.

I've also noticed that all the basemaps are inaccessible through the basemap widget (just shows as loading)... not sure if this is a related issue or not.

What we've tried so far:

white-listed ESRI-specified URLs thinking that it might allow accessibility to these features
re-indexed the portal

Has anyone experienced this and know how to fix it?


texts: Calculating Degrees Minutes Seconds (DMS) to nearest tenth of second in ArcGIS for Desktop?

texts:
Is there a simple way to calculate coordinates in DMS to the nearest tenth second?  I need DD MM SS.s (not DD MM SS.sss).  By default, ArcMap calculates out to 3 decimal places.  I just need it rounded to the 1st decimal place.

I have ArcMap 10.2 and XTools Pro.


texts: Getting closest point from table in PostGIS?

texts:
I have a table of pgrouting nodes, structured like this:

CREATE TABLE ways_vertices_pgr
(
  id bigserial NOT NULL,
  cnt integer,
  chk integer,
  ein integer,
  eout integer,
  the_geom geometry(Point,4326),
  CONSTRAINT ways_vertices_pgr_pkey PRIMARY KEY (id)
)

PGrouting requires the node id, before it can route in-between two locations. My Input are two points(i.e. with Lat & long).

How do I get the node id from the Lat-long? How do I select the closest node?


texts: Converting ESRI SceneLayer Package (.slpk) to Wavefront .OBJ, or other common 3D modelling software format (FBX, STL, etc)?

texts:
ArcGIS.com has some great Creative Commons data I'd like to work with, but the .slpk format it's available in is not well supported in the tools I use (QGIS and blender, and I can't really afford an ArcGIS licence).

How can I convert this 'ESRI SceneLayer Package' (.slpk) JSON vertex/material and texture data into a more common format, like Wavefront .OBJ, .STL or .FBX?


texts:
I made a plugin what has multiple combobox but i have problem with the information load from the second combobox it depends on the combobox 1st.
In the first combo box, i can choose layers from the loaded ones.

    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
         layer_list.append(layer.name())

    self.dlg.reteg.addItems(layer_list)

And the problem is that, in the second combobox i want a list from the selected layer's attribute names.


texts: ESRI Flex Viewer 2.4 - How Can I Remove The Header (or UI Component that holds the HeaderController Widget)

texts:
I am trying to make some modifications to the general layout of the Flex Viewer.  As part of these modifications, I am trying to remove the Header altogether.

I removed the HeaderControllerWidget from the config.xml - thinking the map should now occupy the entire browser.   But now I just get empty white space where the HeaderController was, leading me to believe some UI component must be defined elsewhere to hold the HeaderController widget.

I have searched through the code trying to find this component with no luck.  Does anyone have any suggestions as to how I can remove Header altogether?  I want the map to occupy the entire browser.


texts: coordinate definitions from nyc.gov open data SpatialLocation field

texts:
The NYC Open Data website has a dataset of street events. The last column is called SpatialLocation, and is defined as a text field. An example is below. It is clearly two points on a physical map. Are these points in some standard axis or is it specific to their data? I'd like to convert them to lat, lng coordinates.

GEOMETRYCOLLECTION (LINESTRING (988207 219050, 986291 215593))


texts:
Is it possible to get the SDE version of layers currently used in a MXD through python?

I can get a list (using arcpy.ListVersions) of all versions in the SDE connection file but I only want the version used in the current mxd. I'm using ArcGIS 10.2.1

import arcpy
mxd = arcpy.mapping.MapDocument(CURRENT)

database = "Database Connections/SDECONNECTION.sde"
versions = arcpy.ListVersions(database)

#Print the versions available to the user
for version in versions:
    print(version)


texts:
Methodology from basic Basic Level for the whole process.
I need a list of zip codes. This is the map: https://eligibility.sc.egov.usda.gov/eligibility/welcomeAction.do?pageAction=sfp
If we drill in, we can see on the map that some areas are shaded and some are not. I would like to now any zip code that touches unshaded areas and any zipcode that touches shaded areas. I want a list for all the biege areas and then a list for all non-beige areas.
Any helping material or tutorial.


texts: Several symbols in one line behind a label?

texts:
Is that possible in qgis 2.0:

If a city have several POIs which I can not draw at their origin location because of the scale of the map or because there are just too much symbols, can I draw them in one line behind, or below a label?


texts: How to create geojson Polygons from Lines?

texts:
Given a series of overlapping lines (MultiLineString) in Geojson how can I create polygons for those portions that form an enclosed a space?

Sample data - https://gist.github.com/maphew/7bb77059574ea714977edae41e4b114c

Operational environment is web browser and/or a server with Node.js. Desktop tools such as ArcGIS and QGIS are not available. Connection to a remote service API is acceptable but prefer to avoid. (Following principle of minimizing external dependencies.)

This is my first significant foray into javascript-land, some hand-holding may be needed. ;-)

What hasn't worked

Seeing OGR DXF Closed Polyline Conversion I tried to filter the geojson through DXF to shapefile polygons using command line ogr2ogr (v2.2.0). Thinking if I could get it to work there I could replicate in node-gdal. A polygon shapefile is created but it has no features. -skipfailures didn't help.


texts: How can I extract street information from an osm.pbf file to Excel?

texts:
I'm working on a project where I have to get all the streets of the city where I live, including colony, suburb, postal code, but I am a little confused about what process do I have to follow.
I have download three files from geofabrik.de (osm.pbf, osm.bz2 and shp) that belongs to my country Mexico.
I have other shapefiles of my city (I opened them in the site https://mapshaper.org/ and only could see the boundaries of my city).
So I would like to know how can I extract the full address including the street, postal code, etc.
I work with Windows.
This is what I have done without results (in my point of view)
I load country.osm.pbf file to OSGeo4w QGIS3 (I am able to see the map), but I do not know how to export the data that I need.
I click on the layers and export to csv, but the generated files are to big and includes a lot of points and information that I do not need.
I also use osm2mssql tool to migrate country.osm.pbf to SQL Server but the generated database contains some points but not all the streets.
I also used osmconvert to transform the pbf file to o5m
Then I run this
.\osmfilter.exe mexico.o5m  --keep="addr:country= and addr:city= and addr:postcode= and addr:street= and addr:state="  --ignore-dependencies | .\osmconvert64-0.8.8p.exe - --csv="@oname @id @lon @lat addr:country addr:city addr:postcode addr:street addr:state" --csv-headline -o=filteredfile.csv

But it only generates about 1000 rows.
I do not know if the geofabrik files includes all the streets of my country or only some points about malls, banks and something like that.


texts: QGIS batch clip does not work. One by one it works fine

texts:
When using a batch clip QGIS gets stuck at some point giving me an unknown error, while just clipping one by one works just fine. Any idea what could be causing this?

EDIT:

log says:

Uncaught error while executing algorithm
        Traceback (most recent call last):
          File "E:/usbgis/apps/qgis2/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 203, in execute
            self.processAlgorithm(progress)
          File "E:/usbgis/apps/qgis2/apps/qgis/./python/plugins\processing\algs\qgis\Clip.py", line 113, in processAlgorithm
            if not engine.intersects(in_feat.constGeometry().geometry()):
        Exception: unknown


texts: How do I add and view *.ecw raster images in QGIS 2.2.0 on Ubuntu 14.04 LTS

texts: How do I make ECW files work in Quantum GIS on my Ubuntu 14.04?


texts: Getting coordinates from geometry in PostGIS?

texts:
This line:

SELECT geom FROM myLayer.geom

will return a hex value.

How would I get decimal values, which is a human-readable coordination?


texts:
How can I use the selected value from a SELECT statement in a PostGIS function.

e.g.
using the geom_way from

SELECT geom_way FROM osm_2po_4pgr ORDER BY ST_Distance(ST_GeomFromText('POINT(21.323329 39.948928)', 4326), geom_way) ASC LIMIT 1;

in

SELECT  ST_Line_Interpolate_Point(geom_way, ST_GeomFromText('POINT(21.323329 39.948928)'));

EDIT: I know this has errors, but what I really want to know is how I could refactor a super long statement like this:

SELECT ST_AsText(ST_Line_Interpolate_Point((SELECT geom_way FROM osm_2po_4pgr ORDER BY ST_Distance(ST_GeomFromText('POINT(39.948928 21.323329)', 4326), geom_way) ASC LIMIT 1), ST_Line_Locate_Point((SELECT geom_way FROM osm_2po_4pgr ORDER BY ST_Distance(ST_GeomFromText('POINT(39.948928 21.323329)', 4326), geom_way) ASC LIMIT 1), ST_GeomFromText('POINT(21.323329 39.948928)',4326))));

To be manageable enough to fit into a Java PreparedStatement


texts: Does anyone know how to make new CRS to reproject some old plans?

texts:
I need to reproject some old Cadastral plans, but I can't find any appropriate CRS on web.

They are in old Hungarian System (Budapest tangent plane stereographic projection) which used in late XIX and beginning XX century for topographical survey.
The central point was marked inside old observatory on Gellert Hill,
Budapest.
The ellipsoidal coordinates of central point: φ=47°29'9.6380"north and λ=19°02'56.9441" east of Greenwich.
Used ellipsoid was Bessel and axis were oriented to south-west.
There's no falsing values (some points have negative coordinates).
I have few points in this stereo-sys, and also in Gauss-Kruger 6 zone
MGI for Balkan projection:

pt.num|_______ Ygk|________Xgk|_____Ystereo |   Xstereo|
  1494| 6584836,63| 5059846,65|     -2826,78| 201119,03|
  1033| 6580079,58| 5057131,79|      1967,28| 203771,57|
  960 | 6582258,96| 5061971,45|      -276,39| 198960,07|

Note that axis are set to Y← and X.↓

Has anyone any hints to solve this?


texts: Error message when joining two dataframes with sf -Error: y should be a data.frame; for spatial joins, use st_join

texts:
In R, I'm trying to count the number of points which fall within a buffer around polygons and add the count result to the original shapefile from which the polygon file was created. I've followed the answer provided by @Guzmán from this question: Counting number of points in polygon using R?. Everything works fine until it comes to adding the result to the original shapefile, which I want to do based on the attribute "ID", i.e. not a spatial join. I get the error message "Error: y should be a data.frame; for spatial joins, use st_join" even though y is a dataframe (at least it is according to is.data.frame).

library(sf)
library(dplyr)

OrigPolys <- st_read("OrigPolys.shp") #Load polygons
Points <- st_read("Points.shp") #Load points
NewPolys <- subset(OrigPolys, X1 != 'NA') #Subset required polygons
Buffer <- st_buffer(NewPolys, 5000) #Create buffer NewPolys
inter <- st_intersection(Buffer, Points) #Find points in NewPolys
int_count <- inter %>%
  group_by(ID) %>%
  count() #Count number of Points in Buffer
as.data.frame(int_count) #Convert to df (not sure if this is required as is.data.frame = T whether included or not)
OrigPolys_Pts <- left_join(OrigPolys, int_count, by = "ID") #Join count to OrigPolys based on "ID" attribute```


texts:
I have a shp file with a grid (polygons). For each cell I have the attribute "Join Count" generating during a spatial joints (It is the number of points in each cells).

I would like to export this shp file as a raster.I tried "feature to raster" but no result so far... With feature to raster tool, I set it like this:

But then as result, the raster cells are not overlapping with my grid:

In blue is the new raster generated.

Am I missing something here?


texts: Extract points/pixels from raster based on % of intersection using PostGIS and PostgreSQL

texts:
I need to overlay a geometry (Shapefile of some city in UK) over a raster (grids of entire UK Map) and extract cells/pixels from raster that intersect with points of geometry.

Current Implementation:
I am able to extract the points using below functions:

ST_PixelAsPoints(ST_MapAlgebra(rast,1,ST_AsRaster(geom,rast,'8BUI',1,0,True),1,'[rast.val]')))

Problem with above, is that it extracts all the points, including those points which just touch the raster during overlay.

Expected Implementation:

I need to extract only those points/pixels which overlap above a certain % on the raster. i.e. I need to extract points based on % of intersection.

In essence , I need to replicate the functionality provided by extract function of R . Also more importantly , I need to replicate the weights option provided by extract function.
https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/extract


texts: Prediction of Sentinel images by SVM classifier

texts:
I'm looking for a way to reduce the time taken by SVM classifier (without changing the dimenisons of my training data or data on which prediction has to be done) to predict the raster image data.
I have used the below classifier to predict the data .
    from sklearn.svm import SVC
    import numpy as np
    clf = SVC(C=4, kernel = 'rbf',gamma= 0.25)
    clf.fit(X, y)

    bands_data = np.dstack(bands_data) #bands_data is the stack of raster images(there are around 11 files) data in a ndarray.
    bands_data =np.nan_to_num(bands_data)
    rows, cols, n_bands= bands_data.shape
    n_samples = rows * cols
    flat_pixels = bands_data.reshape((n_samples,n_bands))
    r = clf.predict(flat_pixels)


texts:
I am experimenting with TileStrata and am new to it.  I have a simple tile server set up just to test as follows:

var tilestrata = require('tilestrata');
var disk = require('tilestrata-disk');
var mapnik = require('tilestrata-mapnik');
var strata = tilestrata();

// define layers
strata.layer('heightmap')
    .route('tile.png')
        .use(disk.cache({dir: './tilestrataCache/heightmap'}))
        .use(mapnik({
            pathname: './bin/heightmap.xml',
            tileSize: 256,
            scale: 1
        }));

module.exports = strata;

My mapnik xml is as follows:

<Map srs="+init=epsg:3857">
    <Style name="heightmap">
        <Rule>
            <RasterSymbolizer/>
        </Rule>
    </Style>

    <Layer name="heightmap" status="on" srs="+init=epsg:4326">
        <StyleName>heightmap</StyleName>
        <Datasource>
            <Parameter name="file">E:\GIS_DATA\SRTM\srtm_30m\srtmvrt.vrt</Parameter>
            <Parameter name="type">gdal</Parameter>
        </Datasource>
    </Layer>

</Map>

This all works and a tile is produced without error, however it is always an empty four band (rgba) file, regardless of where I try in the data set.  I have tried 'raw' SRS parameters such as "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over" but I get the same result.  I am hoping to get a grayscale single band 16-bit raster out of this, which should be possible with the mapnik imageType setting (though as yet I don't know where/how to set that as I've fallen at the first hurdle!).

Any suggestions as to what I am missing or doing wrong to get a tile displaying properly (grayscale 16 bit or otherwise!).

EDIT:
I have just tried the Mapnik Hello World example (updating the XML to suit Mapnik v3.x and adding a new layer in the relevant node file) through my Tilestrata stack and get the same result.


texts: Showing some feature in text label in Print Composer from Print Layout itself in QGIS 3.4.5

texts:
I was wondering if there is some expression that I can insert in a text label QGIS print layout that will allow me to print for example coordinates of the map inserted in the print layout itself.

Other features that I'd love to display automatically are:

Before said converted coordinates in another CRS
The scale of the map in print composer


texts: Points displacement - capture new coordinates?

texts:
I have a vector layer of points - say, rooms for rent in a city. Definitely, several rooms are up for rent in one building, so these would share the same coordinates. I need all my rooms to have unique coordinates (say, rooms 1 and 2 in the same building just a few meters apart), so I try points displacement.

Points displacement seems to work visually, but how will I get the new coordinates of the displaced points? "Add geometry column" still exports the original coordinates.

/Edit: I just learned the term, so: what I want is to get the coordinates of the digitized points. I end up with a ring of points (rooms) around the original point (building), but the coordinates of those rooms still end up the same as the original building.


texts:
I would like to plot a raster in WGS84 projection using R based on its attribute value.

Here is the whole procedure

library(maptools)
library(rgdal)
library(raster)

shp <- readOGR("final_result.shp")
cell_size <- 0.012
lon_min <- 10.8; lon_max <- 18.2; lat_min <- 42; lat_max <- 48
ncols <- ((lon_max - lon_min)/cell_size)+1; nrows <- ((lat_max -  lat_min)/cell_size)+1
snap1 <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=0.012, crs="+proj=longlat +datum=WGS84")
snap2 <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=0.012, crs="+proj=longlat +datum=WGS84")
snap1v <- shp[shp$snap == 1,]
print(snap1v)
     class       : SpatialPolygonsDataFrame
     features    : 1
     extent      : 10.9809, 18.18933, 42.07696, 47.89763  (xmin, xmax, ymin, ymax)
     coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
     variables   : 8
     names       : snap,CO,NH3,NOx,SOx,PM25,PM10,VOC
     min values  :  1, 0.0510012180001,0,1.25426065201,1.671099876,0.049, 0.0735,0

snap2v <- shp[shp$snap == 2,]
r_snap1 <- rasterize(snap1v, snap1, fields=c("CO", "NH3", "NOx", "SOx", "PM25", "PM10", "VOC"))

And description of the raster file:

r_snap1
   class       : RasterLayer
   dimensions  : 500, 617, 308500 (nrow, ncol, ncell)
   resolution  : 0.012, 0.012  (x, y)
   extent      : 10.9809, 18.18933, 42.07696, 47.89763   (xmin, xmax, ymin, ymax)
   coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
   data source : in memory
   names       : layer
   values      : 1, 1  (min, max)
   attributes  :
       ID snap         CO NH3      NOx    SOx  PM25   PM10 VOC
        1    1 0.05100122   0 1.254261 1.6711 0.049 0.0735   0

How can I plot the attribute value e.g. CO?
Just using plot(r_snap1) plots value 1 (ID).


texts:
I know in OpenLayers 2 there is the concept of "base layer", which can be set by map.setBaseLayer(layer). But I am not sure which effect this has and if there is an equivalent in OpenLayers 3.

Can anyone help me?


texts: Detect whether there is a spatial polygon in a spatial extent

texts:
Given a spatial extent (for example cbind(c(0,0),c(10,10)), I want to see whether there is a part of a SpatialPolygons falls in this extent?

How can I do this?


texts: generate vector lines of planting from a raster uav image

texts:
I have a high quality drone image and I want to generate the planting line from raster. I will transform rgb to false ndvi or others index. The vector lines must  go through the green areas of the raster (picture)

Is there a plugin or a combination of them to allows me to do it?


texts: QGIS - how to see the minimum and maximum elevations in contour layer

texts:
I am using QGIS 2.12 on a Macbook Pro running OS X Yosemite 10.10.5.

I have a contour layer representing elevations. How can I find what the maximum and minimum elevation values in the layer?


texts: Sentinel API query, rows exceeds maximum number

texts:
I am trying to query Sentinel1A data using the Sentinel API, by using the following scripts.

from sentinelsat.sentinel import SentinelAPI
api = SentinelAPI('username', 'password',
                          'https://scihub.copernicus.eu/dhus/')
api.query(location_str, '20150101', '20150110', producttype='GRD')

I am getting the following error

  sentinelsat.sentinel.SentinelAPIError: (HTTP status: 500, code: None) # HTTP Status 500 - Request processing failed; nested exception is java.lang.IllegalArgumentException: parameter rows exceeds the maximum value (100)

Scihub and API Hub recently updated to limit the number of rows to be returned in a single query to 100, as posted here might be the reason for this error.

How to make OpenSearch API query to return less than 100 rows ?


texts: How to combine sfc objects from R package sf

texts:
Using the R package sf, how does one combine sfc objects? For example, given the following code, how would one create a single sfc object sfc12 that includes the geometries from both sfc1 and sfc2? (length(sfc12) should be 2.)

library(sf)
pt1 = st_point(c(0,1))
pt2 = st_point(c(1,1))
sfc1 = st_sfc(pt1) # An sfc object
sfc2 = st_sfc(pt2) # Another sfc object
# sfc12 = ?

Some approaches that don't work:

sf_sfc(sfc1, sfc2)
# Error in vapply(lst, class, rep("", 3)) : values must be length 3,
# but FUN(X[[1]]) result is length 2

sfc1 + sfc2 # Seems to add the points coordinate-wise.
# Geometry set for 1 feature
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 1 ymin: 2 xmax: 1 ymax: 2
# epsg (SRID):    NA
# proj4string:    NA
# POINT(1 2)

rbind(sfc1, sfc2)
# [,1]
# sfc1 Numeric,2
# sfc2 Numeric,2


texts: Split a field in two separate fields in QGIS

texts: I have a huge sewer network with pipe information, labelled for example as "pvc400" (where material=Pvc and diameter=400). Is it possible to split up the information into text and numbers only, so that i create two new columns, 1 with information about material "pvc" and one with the diameter "400"?


texts: ExecuteSQL on GDAL dataset fails when ORDER BY clause contains more than 1 field

texts:
using GDAL 3.0.4 on windows.

Running an ExecuteSQL command on various datasets (shapefiles, geodatabases, ...) works in almost all cases except when i have an SQL query with an ORDER BY containing multiple fields.

For example either of these succeed:

"SELECT * FROM cities ORDER BY CITY_NAME"

"SELECT * FROM cities ORDER BY COUNTRY_NAME"

But this fails:

"SELECT * FROM cities ORDER BY CITY_NAME, COUNTRY_NAME"

The problem in the failure case, is no error is generated, and the returned OGRLayer is not null.  The OGRLayer is just in a bad state, and calling any methods on the OGRLAyer causes a crash. (eg. null pointer access)    For example, calling GetFeatureCount() on the returned OGRLayer causes crash.

It is not clear from the GDAL documentation if multiple fields in ORDER BY is supported or not.  But even if not, it should return an error and a null OGRLayer.


texts: Error when opening attribute table of a feature layer with a programmatically generated join

texts:
I'm creating a query table with IQueryDef and then using MemoryRelationshipClassFactory to create an inner join between a feature layer and the query table. Both the query table and the join are created without incident and work correctly in ArcMap, (only the features with matching records in the query table are displayed on the map).  However, when I attempt to open the attribute table for the feature layer in ArcMap I get the following error:

Could not load data from the data source.  If you can correct the problem, press the refresh button to reload data.  Possible problems can include bad network connection, invalid field, etc.
Field is not editable.

The table then opens with empty rows.  I'm not sure what the "Field is not editable" message is telling me.  I thought perhaps there was something wrong with the query table, but I can create the table and add it to ArcMap as a stand alone table and it opens correctly. It seems there is something in the way I am doing the join that is causing trouble with the feature layers Display Table.

Below is the relevant code:

    Dim qd As IQueryDef = ftrWksp.CreateQueryDef
    qd.Tables = "swCBCleaning"
    qd.SubFields = "*"
    qd.WhereClause = pre & "CompletedDate" & suf & " >= date '2012-11-08 00:00:00'"

    Dim qn2 As IQueryName2 = CType(New TableQueryName, IQueryName2)
    qn2.QueryDef = qd
    qn2.PrimaryKey = "swCBCleaning.OBJECTID"
    qn2.CopyLocally = False
    Dim dsName As IDatasetName = CType(qn2, IDatasetName)
    Dim dataSet As IDataset = CType(ds.Workspace, IDataset)
    dsName.WorkspaceName = CType(dataSet.FullName, IWorkspaceName)
    dsName.Name = "swCBCleaning"
    Dim name As IName = CType(qn2, IName)
    Dim tblJoin As ITable = CType(name.Open, ITable)

    Dim memRelClassFactoryType As Type = Type.GetTypeFromProgID("esriGeodatabase.MemoryRelationshipClassFactory")
    Dim memRelClassFactory As IMemoryRelationshipClassFactory = CType(Activator.CreateInstance(memRelClassFactoryType), IMemoryRelationshipClassFactory)

    Dim relClass As IRelationshipClass
    Dim dispRelClass As IDisplayRelationshipClass

    relClass = memRelClassFactory.Open("Join", geoFtrLayer.DisplayFeatureClass, "GlobalID", CType(tblJoin, IObjectClass), "FeatureGlobalID", "forward", "backward", esriRelCardinality.esriRelCardinalityOneToMany)
    dispRelClass = CType(geoFtrLayer, IDisplayRelationshipClass)
    dispRelClass.DisplayRelationshipClass(relClass, esriJoinType.esriLeftInnerJoin)


texts: Adding a raster layer on map canvas in PyQGIS

texts:
canvas = QgsMapCanvas()
canvas.show()
canvas.setCanvasColor(Qt.white)
canvas.enableAntiAliasing(True)
vlayer = QgsVectorLayer('testdata/airports.shp', "Airports layer", "ogr")
canvas.setExtent(vlayer.extent())
canvas.setLayers([vlayer])

tms = 'type=xyz&url=https://tile.openstreetmap.org/{z}/{x}/{y}.png&zmax=19&zmin=0'
layer = QgsRasterLayer(tms,'OSM', 'wms')
canvas.setLayers([layer ])

Above code works fine with vector layer but when I am trying to QgsRasterLayer
show cyan colour.

How can I add raster layer on map canvas?


texts: Generating random points from classified raster layer?

texts:
I need to generate random points using a raster layer, which is classified into 5 classes of TPI (terrain position index). This should be a stratified random sampling design, with an equal amount of points on each class and ideally a minimum distance between the points.

I prefer to work in ArcMap. Does a tool for this exist? I heard Hawths tools could work, but I have v.10.3 so maybe not.

You can see one of the sites (classified TPI) here:


texts: updating Z coordinates of points in polylines

texts:
I have a polyline feature dataset with Z coordinates. I'd like to change the Z values (e.g.) by adding 10 meters.

I am trying the following code, but no changes happen. It seems the code is able to enter each part of each polyline and print the Z values plus 10, but I don't know what I have to do to update the Z coordinates.

Here is the code:

import arcpy

gdb = r'C:\Umberto\network_analyst\umbeNETWORK\OfficeBuilding.gdb'
feature = gdb + r'\Transportation\FloorLines'
fields = ["Floor_Number", "SHAPE@"]

# start edit session
edit = arcpy.da.Editor(gdb)
edit.startEditing()
edit.startOperation()

# Create field name with the proper delimiters
#
whereclause = """%s = 3""" % arcpy.AddFieldDelimiters(feature, fields[0])

with arcpy.da.UpdateCursor(feature,fields, whereclause) as cursor:
    for row in cursor:
        polyline = row[1]
        for part in polyline:
            for pnt in part:
                pnt.Z = pnt.Z + 10
                print pnt.Z
        cursor.updateRow(row)

# stop esit session
edit.stopOperation()
edit.stopEditing(True)


texts: Running ArcGIS model at specified time or with delay?

texts:
I am looking to start running an ArcGIS model built in ModelBuilder in the future.

Start at a specific time.

or

Start in x hours.

or something similar.

When I leave the office I may have data processing in GDAL that I know will run for 3-6 hours.

I then want to process the outputs in ModelBuilder but I need to set it to start in say 8 hours.


texts: Inputting bookmarks from TXT file using PyQGIS

texts:
I'm trying to import spatial bookmarks from the QGIS Python console, adapted from this thread and this blog, for the text file I created of exported bookmarks from Arc
import sqlite3
dbPath = QgsApplication.qgisUserDatabaseFilePath()
db = sqlite3.connect(dbPath)
cursor = db.cursor()

srcFile = "D:\\Folder\\filename.txt"

with open(srcFile,'r') as infile:
    bkm = infile.readlines()

for b in bkm:
    name,project_name,xmin,ymin,xmax,ymax,EPSG = b.split(',')
    cursor.execute("""INSERT INTO tbl_bookmarks(bookmark_id, name, group, xmin, ymin, xmax, ymax,projection_srid)VALUES(NULL, str(name), str(project_name),float(xmin),float(ymin),float(xmax),float(ymax),int(EPSG))""")

But I get the following error:

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "", line 1, in
File "", line 13, in
sqlite3.OperationalError: near "group": syntax error

Which I'm not sure what to make of.
I've tried various iterations of name capitalization (e.g. "Group" instead of "group") consistent with the column names

But with no joy.
Any thoughts?
PS. My other thought is to write to XML out of Arc instead of text and import via the QGIS Spatial Bookmarks Manager GUI:

But prefer the console-based approach.


texts: What is proper way to get multiple fields from da.SearchCursor?

texts:
I just buried a couple of hours not realizing that I can't do this:

with arcpy.da.SearchCursor(fp_in, ['OBJECTID', 'name']) as cursor:
  objectid = [row[0] for row in cursor]
  name = [row2[1] for row2 in cursor]
  # now do some stuff with those two lists

the objectid variable comes out as a list with numbers from the feature, as it should.  However the name field comes up [] (empty).  Apparently I don't know with blocks or list comprehensions because it seems I destroyed the cursor object after first use?  I am confused.  name IS a field and it exists.  I solved this by doing the ugly initiate a blank list and append:

objectid = []
name = []
with arcpy.da.SearchCursor(fp_in, ['OBJECTID', 'name']) as cursor:
  for row in cursor
    objectid.append(row[0])
    name.append(row[1])

But there has to be a prettier way right?


texts:
I am trying to get all layers from gdb file.I am using file geodatabase API. I am unable to find any help regarding how to get all layers and then read that layers.

Right Now i am using
 Geodatabase geodatabase = Geodatabase.Open("../../samples/data/Querying.gdb");

//   String allTables = geodatabase.DataSetTypes.;

I have tried several things like get datasetypes and from there get Tables but no luck.


texts: With multi-line labels, how to show the 2nd line only for labels that need it?

texts: Mac OSX 10.6.8 QGIS 1.8.0 : I've two columns in the attribute table, name, and name2 for those labels that need to be 2 lines; multiline labelling is checked. But when I use "name"||'\n'||"name2" I only get labels when there is a printable entry for both name and name2. Can you spot what I'm missing to show the 2nd line only for labels that need it?


texts: Displaying GeoServer WMS GetFeatureInfo in ArcGIS Pro?

texts:
I recently set-up a GeoServer WMS endpoint with GeoWebCache.

I can get the WMS layer to display fine in ArcGis Pro 1.2.

When I click on the map to get feature information then I can see the request in Fiddler being made by ArcGis to the GeoServer.

The response contains the correct HTML but ArcGis never displays it in a pop-up.

This is not an issue in QGis.

Does anyone know why ArcGis Pro is not displaying the contents of the GetFeatureInfo response?


texts: GeoJSON filtering based on properties not working

texts:
I am getting the GeoJSON URL from GeoServer and the data is getting displayed on the map using Leaflet. Now I need to display a few of the features based on the property "status" whether it is pass or rescan.

But I'm getting the following error:

Uncaught TypeError: Cannot read property 'status' of undefined

Below is my code:

var Dupont = new L.GeoJSON.AJAX("http://localhost:8080/geoserver/ABHI/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=ABHI:Dupont_new&maxFeatures=50&outputFormat=application%2Fjson",{dataType:"json"});
Dupont.addTo(map);

var stataupass=L.geoJson(Dupont, {
 filter: function(feature, layer) {
    return (feature.properties.status)=="rescan";
  }

}).addTo(map);

link to Geojson data (fiddle)


texts: Converting CLC local urban atlas vector files to raster

texts:
I am try to convert vector files to raster files from the Corine land use data set which provides vector files in their urban atlas. As per:
https://land.copernicus.eu/local/urban-atlas
For example, I downloaded the files for Wien from here:
https://land.copernicus.eu/local/urban-atlas/urban-atlas-2012?tab=download
In particular, I am interested in rasterizing the field CODE2012, which provides the integer code for the land classification (ITEM2012 and OBJECTID could possibly be an option but the most direct route is CODE2012).
I tried using the QGIS->Raster->Conversion->Rasterize which uses GDAL but I keep getting empty geotiffs or errors. I have also tried manually with gdal_rasterize which is what QGIS uses:
gdal_rasterize -l AT001L3_WIEN_UA2012_Boundary -a Pop2012 -ts 1.0 1.0 -a_nodata 0.0 -te 1843982.5417865077 6166351.096614478 1848539.9284615591 6169446.97652822 -ot Float32 -of GTiff ./AT001L3_WIEN_013/Data/AT001L3_WIEN_013.gpkg output.tif

I have tried many different versions of the above, that is just for an example.
I have looked at vector to raster conversion posts and understand the concept. I need to figure out how to get this particular data set conversion working.
Any help is appreciated.


texts: Adding script tool to ModelBuilder in ArcMap?

texts:
I have create a script and want to add it to the model I have provided. The code is:

import sys, string, os
import arcpy
import numpy
Table = arcpy.GetParameterAsText(0)
Field1 = arcpy.GetParameterAsText(1)
Field2 = arcpy.GetParameterAsText(2)

arr = arcpy.da.TableToNumPyArray(Table, Field1)
arcpy.AddField_management(Table, Field2 , "DOUBLE")
arcpy.CalculateField_management(Table,Field2,(numpy.mean(arr[Field1])))

out = arcpy.SetParameterAsText(4, Table)

However I want to select a field of a feature and calculate the mean value of that field and with the result of the mean continue the rest of the analysis in the model. I have create 3 Parameters, that 1th is a table and the input field and output field. However, when I import the script tool in the model , there is no output feature that I could link it to the rest of the model. How can I do that?


texts:
What is the convention for acknowledging the source of the base map used in a map generated in QGIS via the OpenLayers plug-in?

I'm a teacher of high school geography students who use QGIS to generate maps from data collected in the field, and their work is submitted to an external organisation for assessment, so academic honesty is especially important.

Students are entering data into .csv files and creating their own shapefiles, which I'm fine with... but they're putting that data over a base map - usually one of the OSM maps available via the OpenLayers plugin (QGIS 2.14). I currently have them simply put 'Base map: [base map name]' within their map info to at least acknowledge the source. See example below:


texts: Use Symbology of another shapefile within the same project

texts: In ArcMap, I can go to Properties/Symbology and click on Import to use the symbolization of another shapefile within the same MXD. Is it possible to do the same thing in QGIS? Since I am importing several files that need the same symbolization (symbol, colors, etc) it is very time consuming to do it for every single shapefile.


texts: ArcGIS hillshade create contour artifact

texts: I created a DEM in ArcMap (10.3) using Topo to Raster.  Then, when I tried to create a hillshade using spatial analyst, result look like the picture below.  Any idea on what causing the contour-like artifact?  The profile graph on the left represents the values of the hillshade along the black line, while the graph on the right represents the DEM values at the same location.


texts: What to set spatial reference to after projecting to GCS_WGS_1984

texts:
I have a shapefile containing some sort of vectors or polygonal chains. This file is in the GCS_CH1903 geographic coordinate system (it is not specified though... if I import it into ArcMap it says Unknown). I want to project it to the GCS_WGS_1984 coordinate system and display it on a basemap(also in GCS_WGS_1984).
So I used Define Projection from the ArcToolbox on my imported layer and set it to GCS_WGS_1984. I think the projection is working, because the size of the lines seems correct, but there is no spatial reference so it was displayed somewhere west of Equatorial Guinea in the sea. I tried Create Spatial Reference but nothing happens if I select either GCS_WGS_1984 or  GCS_CH1903.

Layer properties:

Projected Coordinate System:    WGS_1984_Web_Mercator_Auxiliary_Sphere
Projection: Mercator_Auxiliary_Sphere
False_Easting:  0.00000000
False_Northing: 0.00000000
Central_Meridian:   0.00000000
Standard_Parallel_1:    0.00000000
Auxiliary_Sphere_Type:  0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree

The information provided in this link is probably necessary but I have can't figure out where to enter it: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/swiss_grid.html


texts:
I'm working with land use rasters in QGIS. For my application I have created rasters with pixels that have a fraction of 0-1 (for 0-100%) for a land use class. So for example in my settlement raster the pixel value is 0,2. The overlapping pixels from the other rasters have e.g. 0,4 for cropland and for 0,4 grassland.
I need to adjust the pixel sums of these rasters to bring them in line with statistical areas.
So I need something, that increases or lowers the values of the Pixels to a point where it is in line with the new sum value I provide, but at the same time avoids values to get smaller than 0 or larger than 1.
My thought was about something, that loops through the raster and increases/decreases the values to a point where the pixel sum is close to my wanted sum.
I have some basic programming skills, but I never really scripted something in Python before.
How can I solve this with a tool or the graphic modeler? Otherwise, I would also go into scripting.


texts: Distance between Paris and San Francisco with PostGis (out of range coordinates?)

texts:
SELECT ST_Distance(gg1, gg2) As spheroid_dist, ST_Distance(gg1, gg2, false) As sphere_dist
FROM (SELECT
    ST_GeographyFromText('SRID=4326;POINT(37 -122)') As gg1,
    ST_GeographyFromText('SRID=4326;POINT(48 2)') As gg2
    ) As foo  ;

--> ERROR:  Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

I understand those points are pretty distant (San Franciso vs Paris, France), but  how can I get the distance calculated between those two?


texts:
I am using 'Esri brochure' to load URL ArcGIS tourism services with the API 'L.esri.DynamicMapLayer'.
I'm looking for a way to know how many layers there are in that service, preferably the layer id.

I.e:

this service:
sampleserver5.arcgisonline.com/arcgis/rest/services/Hurricanes/MapServer

has two layers

sampleserver5.arcgisonline.com/arcgis/rest/services/Hurricanes/MapServer/0
sampleserver5.arcgisonline.com/arcgis/rest/services/Hurricanes/MapServer/1


texts: How to count legend items in QGIS print composer

texts:
Is there a way to count legend items displayed in the print composer ? I would like to adjust by expression the number of columns of my legend based on the number of legend items.
Example : if there are 6 or less legend items I want my legend to be displayed on 3 columns. If the number of items exceeds 6 then I want my legend to be displayed on 4 columns.
You can find below a screenshot for illustration purposes which features 7 legend items with a column count of 4.
Thanks


texts: Making curved line between points by using QGIS Geometry Generator

texts:
I am testing different tree crown drawing options in QGIS. At the moment, I have found this option to be the most suitable for me.
Code:
convex_hull(
            combine(buffer(make_point($x-("west"/2),$y),"west"/2),
            combine(buffer(make_point($x+("east"/2),$y),"east"/2),
            combine(buffer(make_point($x,$y+("north"/2)),"north"/2),
                    buffer(make_point($x,$y-("south"/2)),"south"/2))))
            )

Problem is that uneven crowns look very good, but even crowns look like rounded squares. I would like lines to curve a bit. Right now they are straight.
I understand how generator draws them, but not enough to make it work like I would like. Info from the code author and how it works. Look for Solution 3.


texts:
I am having trouble publishing my US Address - Dual Ranges Address Locator to our ArcGIS Server. I've created it with a very small sample of 7 streets around my local University just to test, with no odd address types (built on US Census Address Features shapefile).

It's properly located (myserver\arcgisserver\arcgisresources\Locators) and works when I reference it in ArcGIS Desktop.

Yet when I publish it to our ArcGIS Server, it gives me an error any time I enter a specific house number. Works finding street intersections and even spits back correct addresses if I reverse geocode. But it will not work searching for a specific address.

Anybody understand why this is happening? I am at a loss. I need this to work for my users.

EDIT: I am working in ArcGIS 10.0 SP 5 and ArcGIS Server 10.01.


texts:
i have installed QGIS with GRASS and QDAL complete...

I followed the instructions for additional QDAL configuration (adding usr/bin to the QDAL path).

When i try to "fetch plugins" my list is short and does not include geocoding.

I have downloaded a plugin from the main repository but it sits on my desktop undetected by QGIS

How do i access these plugins?


texts: Popup shows the data twice when exporting a specific layer to OpenLayers with QGIS2Web

texts:
I've got a problem when exporting a specific layer to OpenLayers with QGIS2Web.

This specific layer is a Parcel or Lot and has subdivision parcels or subdivision lots within. The main parcel has data as area, owner, location... and same data for each subdivision lot.

However, when I export it to OpenLayers and click on a parcel it shows me the Main parcel data and the subdivision lot data, so because they are the same owner and same location and I see the very much the same information except for the area.

But, I export it to Leaflet it works as I want it to work. When I click on a parcel it shows me just one set of data. The subdivision lot.

I really want to use OpenLayers so I would like to know if someone has a solution for this.


texts: How to cut a 'tiled image' from a raster dataset using GDAL?

texts:
I am working on an app that needs to create tiles (jpeg) from a given raster dataset.  My initial inclination was to use gdal2tiles.py, but its performance seems like it could be improved.

That leads me to ask the following question:  What would be the fastest way to cut a 'tiled image' from a raster dataset using GDAL?  And by 'tiled image', I mean just a simple jpeg or png.

In my prototype, I use the MapServer C# MapScript bindings to do the job.  That is, I create a map object, loop through all the tile bounds I need to cut, set the map objects extent, and then save the resulting image.  The performance of this approach is significantly better than that of gdal2tiles.py, but I am wondering if I use GDAL directly if I can't make it even faster.  Can anyone suggest a similar workflow with GDAL methods?

EDIT:  After some more research today, I found the answer right in front of me.  If you have downloaded FWTools, the csharp\apps folders contains a number of classes to demonstrate the C# GDAL bindings.  In my case the GDALRead.cs and GDALReadDirect.cs were what I was looking for.


texts:
I try to create a multiple lines label for this this polygon:

I red

How to create multiline labels in QGIS?

QGIS, Multiple lines and wrap on character?

and tried to use the "wordwrap" function with this expression:

wordwrap('high sensivtity area - red',10  )

and in the output preview (down the Expression string builder) it look fine and this is the result i'm seeking - but when i choose OK - nothing change in the map.

I don't want to use the option to put a space in the 'Wrap on character' section, because i don't want each word i separate line


texts: Selecting more than one country to be highlighted

texts:
I can get a one country within a world shapefile to highlight, but when it comes to adding additional countries this just turns all countries back on

I am using query builder to achieve highlighting Mali, but really want to highlight Burkina Faso and Niger

If I do "admin" = 'Mali' this works but if I do "admin" = 'mali' AND 'burkina faso' this doesn't work

I have also tried doing

"admin" = 'mali'

"admin" = 'burkina faso'

This doesn't work either


texts: Calculating proportion of area that is land within bounding box using Python

texts:
I am trying to find the proportion of land in an area specified within a bounding box (in python). The bounding box is given by min/max latitude & longitude coordinates.
The initial solution I came up with was to sample N random points within the bounding box and check (for each point) whether it lies within a landmass or not. Specifically, I checked the point within a landmass shapefile (from here: https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-land/). It doesn't give an exact area answer but an estimate.
This solution works but it is really slow (even with multithreading) so I was wondering if there is a better way to do this? I was thinking that you could isolate a landmass polygon within the bounding box from the shapefile and calculate the area within that. But I am very new to the osgeo package and shapefiles in general so I am not sure where to start with this problem.
Is there another way that I could do this?
Here is my reference code and some pictures of my initial solution (some code is missing like the bounding box calculation code, I can add it at request if needed):
import numpy as np
from osgeo import ogr
from itertools import repeat
import multiprocessing as mp
import matplotlib.pyplot as plt

def is_pt_land(lat, lon):
    shape = ogr.Open('./ne_10m_land_shape/ne_10m_land.shp', 0)
    shape_layer = shape.GetLayer()
    geo_ref = shape_layer.GetSpatialRef()
    point_ref = ogr.osr.SpatialReference()
    point_ref.ImportFromEPSG(4326)
    ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)
    #Transform incoming longitude/latitude to the shapefile's projection
    [lon,lat,z]=ctran.TransformPoint(lon,lat)

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

    # check if point exists within shape
    shape_layer.SetSpatialFilter(pt)
    return len(shape_layer)>0

def is_land_thresh(bounds, iters=5000, thresh=0.80):
    xpts = np.random.random_sample((iters,))
    ypts = np.random.random_sample((iters,))
    xpts = (abs(bounds[0][0]-bounds[1][0])*xpts) + bounds[0][0]
    ypts = (abs(bounds[0][1]-bounds[1][1])*ypts) + bounds[0][1]

    with mp.Pool() as pool: # applying the is land function to all the random points
        truth_arr = pool.starmap(is_pt_land, zip(xpts, ypts))

    truth_arr = np.asarray(truth_arr) # creating a visual of the truth array to generate an estimate
    plt.figure()
    plt.scatter(ypts[truth_arr], xpts[truth_arr], c='green', s=2)
    plt.scatter(ypts[~truth_arr], xpts[~truth_arr], c='blue', s=2)
    plt.show()
    return (np.sum(truth_arr)/iters) >= thresh

# Testing bbox and land threshold calc
loc = GeoLocation.from_degrees(26.062951, -80.238853) # create center point
bbox = loc.get_bbox(50)   # get bounding box with 50km from center point (100km x 100km roi)
is_land_thresh(bbox)      # determine if area is mostly landmass

Here is a photo of the scatter plot calculation:

Here is the landmass bounding box denoted by the min lat/long as the SW coord., the max lat/long as the NE coord., and the center coordinate:


texts: Why does ArcGIS ModelBuilder tool Iterate less than expected Number of Times?

texts:
I have an ArcGIS ModelBuilder tool (with an iterator in the workflow) and the tool iterates through approximately 30 features then stops.  I have checked that:

It is not instructed to iterate only a certain number of times; and
That the input feature class does not have a selection on it.

Has anyone else experienced this behavior?

Additional details :

The iterator churns through a dataset and creates a new FC by unique ID.  The next steps are complex and include building terrains.  The behavior is usually corrected when I shut down ArcMap and open up a new MXD.


texts: Displaying specific country into map using countries.geojson file?

texts:
I was trying different ways to display one country using countries.geojson file, but I unable to display.
I tried like this, but map is not displaying....

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

var vector = new ol.layer.Vector({
    source: new ol.source.Vector({
        format: new ol.format.GeoJSON(),
        url: 'https://raw.githubusercontent.com/openlayers/ol3/6838fdd4c94fe80f1a3c98ca92f84cf1454e232a/examples/data/geojson/countries.geojson'
    }),
    style: function(feature, res){
        if(feature.get("name") == "Germany"){
            return new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'red',
                    width: 2
                })
            });
        }
    }
});


texts: How to use coordinates search bar in QGIS?

texts: When using coordinates search bar in QGIS, it sends me to a completely different location. The GPS coordinates style is 37.381096, 14.561191 and the project is set to EPSG: 3857 with OpenLayers basemap. What am I doing wrong?


texts: A GeoTIFF heightmap from SRTM embeds vertical planes along lat/lon grids in its 3D representation

texts: I've downloaded an SRTM GeoTIFF from Open Topography and converted it into the grayscale image before I imported it into Blender to improve the quality. It is apparent that vertical planes are cutting through the heightmap at places where I guess the lat/lon grid runs. It must have been there from the first in the GeoTIFF and probably be stood out through the enhancement procedure. How can I remove them without having to paint the lines or planes with a brush?


texts:
I'm using QGIS 2.6.1 brighton, and I have a problem in computing the area of each layer.

I was actually working on the map of our municipality and I would like to compute the total area of the residential places etc.

The problem is, it actually just gives me 0 when I compute the areas.

I tried using the prs 92 since I'm from the Philippines, but the base map of Google Earth I'm using has a different projections.

What should I do?


texts: Setting default parameter/style for new layers in QGIS?

texts:
For each new layer I create in QGIS, I have to set the style parameters again.

Can I set the default parameters for each new layer?


texts: Reproject gif.gz raster to WGS84 in QGIS

texts:
I am a steep learning curve into QGIS. I have a raster of global fire maps from the modis site: http://lance-modis.eosdis.nasa.gov/imagery/firemaps/firemap.2013221-2013230.8192x4096.gif.gz

I would like to overlay this raster with a shapefile showing Global administrative boundaries GAB. I have both files in WGS-84. However, I have noticed that the fire raster is offset south-east of the WGS84 GAB shapefile and the two Do not overlay.

Kindly assist


texts: How do I compact a file geodatabase with ArcPy?

texts:
I am preparing some data using FME, which I end up writing to an Esri file geodatabase. Unfortunately, FME doesn't appear to update the written feature classes' spatial indices. This results in display problems when I later inspect the written data with ArcMap (or other ArcGIS desktop tools).

It appears that there are two seemingly easy solutions for this problem:

Recreate any spatial indices after writing data.
Compact the file geodatabase after writing data.

Compacting seems like the easier choice, because it can be done per geodatabase (vs. recreating indices per feature class).

I don't want to compact the geodatabase manually; ideally I'd want to compact the geodatabase inside my FME workbench, using a shutdown Python script. However, there appears to be no ArcGIS geoprocessing tool for compacting a geodatabase (only one for compressing, which is not the same thing).

Does anyone have a solution here?

Updating the spatial indices directly in FME; or
compacting a file geodatabase using ArcPy (in an FME shutdown script).


texts: Count number of pixels satisfying different conditions in an online layer within a polygon

texts:
Raster file-> Terrain: Slope in Degrees in ArcGIS Online, by ESRI

Shapefile-> US Census Tracts in Census 2000, downloaded from NHGIS

I would like to have more information about the distribution of slope within a certain census tract (each of them is a polygon). In particular, how many pixels satisfy different conditions (e.g. slope between 0-5%, between 5-10%, 10-15, >15) and also calculate the range of slope within the polygon.

I've done something like this before using Con and Zonal Statistics as Table but it meant having use the Python console to create a different raster for every possible category and that would take too long (also lots of HDD space). Also, I don't know how to do this with ArcGIS Online Layers since they are huge (this one is about 1.5TB), so maybe there is another option.


texts: calling gdal_merge.py into another python script running GDAL processes

texts:
I want to use gdal_merge.py to merge a series of .tif files prior to clipping the merged file to a shapefile boundary, but want to do so as part of another larger script that will execute a number of other processes.  I am a python/GDAL newbie and am not sure how to go about calling gdal_merge.py into another script.

I cannot merely run gdal_merge.py on its own as it is but one step in a script that will hopefully execute a number of processes. Any thoughts on the best way to do this?


texts: DetectFeatureChanges_management gives out error 000366 "Invalid geometry type"

texts:
I am trying to compare county borders in two shapefiles and create a third showing how (if at all) those borders changed. To do so, I'm using arcpy.DetectFeatureChanges_management.

I have tried the following code:

# Detect changes in county boundaries over time
import arcpy

#Set working directory
cd = r"C:\Users\folder"
arcpy.env.workspace = cd
arcpy.env.overwriteOutput = True

#Set local variables
file1 = r"county_shapefiles\county_2000.shp"
file2 = r"county_shapefiles\county_2005.shp"
output = r"county_shapefiles\county_changes.shp"

#Create changes shapefile
arcpy.DetectFeatureChanges_management(file1, file2, output, "500 Feet")

but I get the error 000366: invalid geometry type. The error occurs twice, so I assume it is related to the two shapefiles I use as input.

I consulted the Help from Esri (http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/detect-feature-changes.htm), but in their example code they use shapefiles (.shp), too, so this should not be the problem.

Any ideas?

Update:
Following ShaunO's suggestion I tried to use the intersect function to do the same. I adapted his code in the following way:

file1 = r"county_shapefiles\county_2000.shp"
file2 = r"county_shapefiles\county_2005.shp"

arcpy.MakeFeatureLayer_management(file1,"lyr1")
arcpy.MakeFeatureLayer_management(file2,"lyr2")

#make an empty set to hold outputs from intersect tool

CountyIDs = set()

with arcpy.da.SearchCursor("lyr1", "COUNTY") as cursor1:
    for row in cursor1:
        with arcpy.da.SearchCursor("lyr2", "COUNTY") as cursor2:
            for rows in cursor2:
                #create sliver with intersect toool for each county
                arcpy.Intersect_analysis(["lyr1", "lyr2"], rows[0]+".shp")
                #Add Interset outputs to set for merge
                CountyIDs.add(rows[0])

#finally, merge all the outputs
arcpy.Merge_management(CountyIDs,"county_changes.shp")

#delete layers to release them from memory
arcpy.Delete_management(lyr1)
arcpy.Delete_management(lyr2)

I have changed some things that produced errors before (no guarantee I changed anything for the better. Still new to ArcPy/Python). Now the code does not produce any error anymore, but it does not finish either. Is this due to a programming error on my behalf, or simply because intersecting close to 4000 counties, and hence adding up 4000 shapefiles, is a very time-consuming or inefficient task?


texts: Choosing coordinate system to use in South Africa?

texts: What are the principal datums and projections used in South Africa including the names of geodetic datums and their projections?


texts: How to convert .OBJ 3D Mesh Georeferenced file generated from OpenDroneMap to .kmz format for Cesium?

texts: I have Open Drone Map Output of a 3D Mesh Object file in .obj extension. I want to convert it into .kmz format to put it into Cesium Ion. How can I convert the files in the picture to a .kmz format file?


texts: Why are colour distorted at large solar zenith angle in satellite image?

texts:
Most Earth observation satellites consider large Solar Zenith Angle as undesirable are therefore do not process them. There are indeed several known image quality issues, such as low signal (inducing low signal-to-noise ratio) and the so called "chromatic artefact" (inducing distorted spectral signatures).

Low signal is due to the cosine relationship between the energy per square meter and the incidence angle.

As far as I understand at this stage, chromatic artefacts are due to the different behaviour of the light going through the atmosphere depending on its wavelength. This includes absorption and diffraction. However, it is not clear to me why it is so difficult to manage by atmospheric correction experts and why this has so much impact on the satellite images (considering that the satellite viewing angle is close to nadir).


texts: Seeking map of communes for Bavaria region before 1930?

texts:
I am looking into data for an article and I am missing a pretty important map. The map relates to this question. I need a map of Bavaria communes, the so called gemeinden divisions, either to buy or freely available, raster or vector format. I would accept any format that can be imported into a GIS software, but it needs to represent the period prior to 1930.

I have spent around 5 hours already and all I could find related to the map I am looking for is in the link to the other question. I mean I found other sites, and old books, but nothing concrete.


texts: Importing CSV file to QGIS from Excel columns with formulas

texts:
I import a CSV file on Qgis.

In excel I see an incription like that ¨54235¨ and in the attribute table in qgis as well as in notepad it looks like that ¨ 54 235¨ *not space but tab.

The column when it was an excel file it was created like text join (e.g: =A11&B12)

The initial columns A & B as well have the same problem (e.g. ¨012¨ apears ¨ 012¨ with a tab.

Do you know a solution to that problem?