texts:
One thing that I can't quite get my head around is Reprojection.

It comes up a lot in questions on Stack Exchange, and there are plenty of How to questions and answers about specific issues using various software.

I have searched the web with "What is reprojection?" and I have had a good read. The ESRI FAQ about Projection Basics is good http://support.esri.com/en/technical-article/000005562 - From this...

The terms 'geographic coordinate system' and 'datum' are used interchangeably, but as noted above, a GCS includes a datum, spheroid, units of measure and a prime meridian. The coordinates for data change depending on the datum and spheroid on which those coordinates are based, even if they are using the same map projection and parameters.

...I feel this starts to explain.

I understand Coordinate Reference Systems. I understand the reality that the earth has an irregular shape and the that different reference systems 'tie' theoretical co-ordinates to real places on the planet surface, and that different CRS are needed around the globe.

This GIS SE post is useful for that Understanding concept of Spatial Reference System?

And this answer with a link to a working manual of projections is a great resource
Learning about Reprojection?

Again this Wikipedia page starts to touch on reprojection https://en.wikipedia.org/wiki/Wikipedia:Graphics_Lab/Resources/QGIS/Reprojection

...to reproject (warp) an existing lat/long project to a specific projection

What is actually happening when you Reproject the data from one projection, to another, projection. Can it be explained in simple language?


texts: QGIS Zonal Statistics for polygons is missing results?

texts:
I am carrying out an analysis using QGIS 2.14.1 where I am dealing with point data that have been randomly displaced by either up to 2 Km or 5 Km. I have generated a buffer of the appropriate radius around each point and now want to obtain mean values for each polygon from a raster layer.

To do this I have used the Zonal Statistics tool but for about 10% of my points no results have been generated even though there is raster data for that location. All the polygons with missing data have a 2 Km radius. The Zonal Statistics tool worked for some of the polygons with a 2 Km radius, just not all. The raster layer I am dealing with has a resolution of 5 X 5 Km.

When I increase the size of the polygon I am able to obtain zonal statistics results, which suggests to me that the size of the polygon may be too small. However, I cannot alter the polygon size because the point can only be within this radius and therefore enlarging the polygon would make the results inaccurate.

I understand that if polygons overlap this can sometimes cause issues. Whilst some of the polygons do overlap, some of them do not. I also tried doing this for one polygon at a time but still did not get results with the ones I tried. I have included an image below where the smaller green polygon is an example of one of the polygons that has missing Zonal Statistics output.

I am aware of multiple other Q&As addressing this type of issue but nothing I have tried has worked so far. The Q&As I have seen are:

Missing buffer values when using zonal statistics in ArcGIS Desktop?
Zonal Statistics as a Table does not calculate statistics for all features in dataset?
Zonal Statistics As A Table missing Features in Results?
Polygons too Small for Zonal Statistics
Zonal statistics results in some NULL values in QGIS?

I do not use Python.


texts: Using CASE to select between two geometry functions?

texts:
I'm trying to cut holes in forest polygons, where there are lakes.

My PostGIS/PostgreSQL code look like this:

CREATE TABLE samp.forest_with_holes AS
  SELECT
    fo.objectid,
    (CASE ST_Intersects(fo.geom, la.geom)
      WHEN TRUE THEN ST_Difference(fo.geom, so.geom)
      ELSE fo.geom
    END CASE;) as geom
  from samp.forest as fo, samp.lake as la;

But I get the error message:

ERROR:  syntax error at or near "CASE"
LINE xx:     END CASE;) as geom

Can anybody help me to see why this is not working - and please suggest a better way :-)


texts: QGIS Group Points Together Based on Multiple Attributes

texts:
I'm using QGIS and I want to group points that are:

within 20m of each other (as I have multiple mobile basestations on the same site)
based on a categorical variable (the mobile operator)

I have 139,742 points that I expect to reduce to around 30,000 by the time I've grouped those within 20m. This 30,000 then should divide between the major mobile operators.

I want to do this without losing all of the data linked with those points. A good way may be to select points based on these conditions and then add in a new column with a shared number to indicate the grouping. Any ideas on how I may do this?


texts:
I am searching for a solution for my current OpenLayers 3 application. We are now in need of adding what seems to be map routing. I have a small area, a building and a parking lot (about 5 or 6 acres) that does not have normal roads populated on it in OSM. I need to be able to add functionality so that a user is able to navigate from their current location (using ol3 geolocation) to another feature on the map.

This seems like it would not be difficult, however there are rooms and hallways in the building that are in this space so I need something that will be able to route lines on a path to help navigate the user from wherever they are. Basically I am looking for something similar to Google maps driving directions but not for driving, or on roads, but inside a building/parking lot. I did look on Google maps, and they have the room outlines, and hallway outlines that I would need. I am not sure how they obtain that data or how I could represent that in my own map and add routing around that.

I looked into OSRM but I am not sure it is the solution I will need.

Does anyone have a book or tutorial to recommend?


texts: Odd behaviour with MULTIPOINT layer in QGIS using delimited files and WKT

texts:
I'm writing some python to create grid lines, either as a series of densified LINESTRINGs or a set of MULTIPOINTs.

The code uses shapely, and I'm writing WKT to a CSV and importing into QGIS(2.18.3) as a delimited file.
I'm running this in a python venv with the following (output from pip freeze)

Fiona==1.7.5
Pillow==4.1.0
Shapely==1.5.17
click==6.7
click-plugins==1.0.3
cligj==0.4.0
munch==2.1.1
numpy==1.12.1
olefile==0.44
psycopg2==2.7.1
pyproj==1.9.5.1
six==1.10.0

The code is as follows:-

from shapely.geometry import MultiPoint, LineString

def graticule(long_degrees=10, lat_degrees=10, segments=1000,
              geometry_type=MultiPoint):
    """
    Generator. Emits a series of geometries
    :param long_degrees: separation of longitude
    :param lat_degrees: seperate of latitude
    :param segments: number of segments (bigger values = smoother)
    :param geometry_type: type of geom (MultiPoint, Linestring)
    :return: yields shapes. Up to client to iterate over these
    """
    # meridians
    for long in range(-180,180,long_degrees):
        coords = []
        longitude = float(long)
        for lat_step in range(0, segments):
            latitude = float(90.0-(lat_step*(180.0/segments)))
            coords.append((longitude, latitude))
        yield geometry_type(coords)
    # latitudes
    for lat in range(-90,90,lat_degrees):
        coords = []
        latitude = float(lat)
        for long_step in range(0, segments):
            longitude = float(-180.0+(long_step*(360.0/segments)))
            coords.append((longitude, latitude))
        yield geometry_type(coords)

def main():
    """
    Makes a 10 degree grid as a WKT CSV file
    :return:
    """
    with open("/tmp/foo.csv","w") as output_file:
        output_file.write("GEOM\n")
        for geom in graticule(10,10,1000,MultiPoint):
            output_file.write("{}\n".format(geom.wkt))
    print("Finished!")

if __name__ == "__main__":
    main()

It works flawlessly when I choose LineString as the geometry_type...

but when I choose Multipoint, I only see things to the east of the Prime Meridian. And the lines of latitude disappear completely.

Does anyone have any clue what's going on here? I've shown using the Bonne projection, but it also happens with plain old wgs84.
I've also tried sorting the multipoints so that they go north to south, but the result is the same :/

No errors were thrown.


texts: Errors fitting a Markvov Random Field GAM with multiple samples per Polygon

texts:
all, after reviewing this excellent tutorial, I'm trying to expand out my use of GAMs for space and time. I'm trying to use the Pennsylvania smoking and cancer data, but am running into a problem when I try and work with multiple values of predictor per county.  Let's consider the data.

library(mgcv)
library(sf)
library(dplyr)
library(SpatialEpi)
library(spdep)

penn_sf <- st_as_sf(pennLC$spatial.polygon) %>%
  mutate(county = unique(pennLC$data$county)) %>%
  left_join(pennLC$data) %>%
  mutate(rate = cases/population*1000) %>%
  st_transform("+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

This is great, but for a MRF, we'll want a neighborhood matrix of the counties - and this data has many rows per county

> penn_sf %>% as_tibble %>% dplyr::select(-geometry)
# A tibble: 1,072 x 8
   county cases population race  gender age        rate    ID
   <fct>  <int>      <int> <fct> <fct>  <fct>     <dbl> <int>
 1 adams      0       1492 o     f      Under.40  0         1
 2 adams      0        365 o     f      40.59     0         1
 3 adams      1         68 o     f      60.69    14.7       1
 4 adams      0         73 o     f      70+       0         1
 5 adams      0      23351 w     f      Under.40  0         1
 6 adams      5      12136 w     f      40.59     0.412     1
 7 adams      5       3609 w     f      60.69     1.39      1
 8 adams     18       5411 w     f      70+       3.33      1
 9 adams      0       1697 o     m      Under.40  0         1
10 adams      0        387 o     m      40.59     0         1
# … with 1,062 more rows

which is fine - I want to look at race, gender, and age all as predictors.  So, let's create a neighborhood matrix for the state, first aggregating at the county level.

#aggregated
penn <- penn_sf %>%
 group_by(county) %>%
  summarize(cases = sum(cases), population = sum(population))

#neighborhoods

nb <- poly2nb(penn,  row.names = as.character(penn$county))
names(nb) <- penn$county

Let's see it!

plot(penn_sf[["geometry"]])
plot(nb, coords =coords, col = "red", add = TRUE)

However, when I then fit the appropriate model...

gam_mrf <- gam(cases ~ s(county, bs = 'mrf',
                         xt = list(nb = nb)) +
                 race + gender + age + offset(population), # define MRF smooth
               data = penn_sf,
               family = poisson)

I get the following error

Error in gam.fit3(x = X, y = y, sp = L %*% lsp + lsp0, Eb = Eb, UrS = UrS,  :
  inner loop 1; can't correct step size
In addition: Warning message:
Step size truncated due to divergence

I've tried a few variants on this based on things I've seen around - using numerics/integers and creating an ID column which is the integer transform of the county factor in penn_sf, a neighborhood matrix created from penn_sf itself, and others - but all yield either this or another error.

1) What's the issue here, and any thoughts on how to fix?

2) Is it correct to use the aggregated polygons, rather than the original one row per data point polygons? It "feels" correct, but that's never a good way to decide things!

BTW: My goal is to then use this model to get fits, sum them up, and see if the observations match total predictions (in addition to looking at it by group) as well as trying a leave-one-out and seeing how well we can predict data for a missing county.


texts: How can I make atlas zoom on different polygons?

texts: I'm trying to create an atlas and my last step is to zoom in on each exported picture over the red polygon. Any suggestions?


texts: Find features that use specific domain in ArcCatalog?

texts:
I have a geodb file with a set of domains.
Now, I need to check which are the features with attributes using a specific domain in arccatalog. Is there any way to accomplish that?


texts: Hide null value attributes in Leaflet popup

texts:
I have marker properties which are fetched from the database and displayed as a table. I want to hide the "reply" feature property when it's null and display it only when there is a text entry against it.

Similarity: hide attribute with null value using leaflet

Although it is a similar question but I am not sure how to get the property labels as shown in the answer.

My code:

    function returnQuery(feature, latlng) {
        return L.marker(latlng).bindPopup(
            "
" +
            "
" +
            "
" + "" + "Query ID: " + "	" + "" + feature.properties.id + "	" + "
" +
            "
" + "" + "Issue Placement Date: " + "	" + "" + feature.properties.q_date + "	" + "
" +
            "
" + "" + "From: " + "	" + "" + feature.properties.q_name + "	" + "
" +
            "
" + "" + "Query/Problem Statement: " + "	" + "" + feature.properties.query + "	" + "
" +
            "
" + "" + "Reply: " + "	" + "" + feature.properties.reply + "	" + "
" +
            "
" +
            "
");
            }


texts: Qgis 3.16.3 - tk package, ModuleNotFoundError: No module named '_tkinter'

texts: I am developing a plugin that uses tk package. When I tried to import tk I got an error: ModuleNotFoundError: No module named '_tkinter'. I tried to reinstall tk but it didn't work. I would say that when i used local python tk package works.


texts: Disable column of Value Table in Python Toolbox?

texts:
Can I disable a column based on another column?

Let's say I have a column to choose a field type (Short, String and so on) and a second column to choose the max length of a string. I would like to disable the second column if I've chosen a short instead of a string.

param = arcpy.Parameter("fields", "Param", "Input", "Value Table")
param.columns = [['String', 'Field Type'], ["Long", "String Length"]]

I'm using ArcGIS 10.3 for Desktop.


texts: How to upload shapefile to PostGIS using cURL and GeoServer REST API

texts:
The GeoServer REST API docs here does not detail how to do this. I'd like to create a PostGIS table and corresponding featureType in GeoServer using the REST API, without having to chain requests.

For now I can upload a shapefile to geosever, but not to PostGIS. I can also create a PostGIS data store in a separate request. Unfortunately, its not possible to use an existing PostGIS data store when uploading a new shapefile. Anyone with experience on how to achieve, please assist.


texts:
QGIS 2.14 worked fine until a couple of days ago (ubuntu 15.10 wily), but now crashes at startup with only this error :

Segmentation fault (core dumped)

sudo qgis doesn't give any error but QGIS don't start.

qgis --noplugins doesn't help.

QGIS was installed from http://qgis.org/debian.

For development, I installed libqwt6. This might have been the source of the issue, but removal of all qwt packages and QGIS packages didn't help.


texts: Purpose of comma in r.mapcalc in GRASS GIS 7.4

texts:
I would like to know the purpose of a comma in r.mapcalc of GRASS GIS. Consider this example -

r.mapcalc "noise = if(landuse96_28m==1 || landuse96_28m==2, roads_buffers, null())"

In the above example I know the statement means if landuse96_28m is equal to 1 OR landuse96_28m is equal to 2 but after this last digit I get confused as to what the  comma(,) stands for as well as the one after the roads_buffers.

If someone would explain to me in a straightforward manner what the comma(,) means I will receive a breakthrough in using r.mapcalc.


texts: Adding id values in attribute table after digitization

texts:
I have digitized land use and land cover polygons and am curious as to why the "id" attribute (presumably automatically assigned and populated) is either NULL or blank for my polygons.

Geometric validation reports feature "id"s, but I cannot relate them to my data. Is there a way to automatically set the "id"s?


texts: saving results of ST_offsetcurve to new table with geometry and fields from original data table

texts:
I am trying to create lines parallel to some road centerlines, using distances from the attribute table specific to each line. I'm trying to figure out how to get it to save to a new table with geometry that I can view in qgis or arcgis. I've imported the source data (which was a shapefile) into a postgres database, and ran this query:

select ST_OffsetCurve(
test_data_.geom, test_data_.buffer_dis, 'join=mitre mitre_limit=5.0')
into test_data_left from test_data_;

It writes out to a new table called test_data_left as specified above, but I can't see this table (or the original for that matter) in either qgis or arcgis, when I try to load directly from the database. I tried to export the test_data_left to shapefile, which it does, but when I try to laod it in qgis it appears to be empty, and arcgis tells me the number or shapes don't match the number of records. I think I am missing a step or two. I'd also like to be able to associate some of the fields in the original table with the ouput but am not sure how to do that in this expression?

UPDATE:

I found this useful thread Creating a spatial PostGIS table as a result of a PostGIS query and am working on adding geometry and ID fields, hoping that will get me part of the way there. Still looking for a way to add some of the original fields and values to the output as well.


texts: I got a featureclass polygon of 4 point polygons that delineate the coverage area of a georeferenced image; these polygons will later be map sheets. Anyway I discovered some of them are not true rectangles and I was wondering if there was a tool or script out there for either Arc or QGIS that could analyse a polygon and reshape it to be a perfect rectangle.


texts: Using onclick event with stream services into ArcGIS Online?

texts:
I have one apartment lat-long.

With this I have 200 meters with some sensor data, and I have published web application with meter data.

I have used stream output with geoevent manager while observing the data only one marker is displayed.

Right now it is showing one real-time data into a single click.

Is it possible to configure with one field means with meter name can be seen with all data in one click?


texts: TypeError: unsupported operand type(s) for -: 'str' and 'int'

texts:
I have 732 temperature rasters (366 min and 366 max temp) that I want to manipulate depending on if the max temperature was greater than a specified and constant crop category temperature.

import os
import arcpy
import sys

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage("Checking out Spatial Analyst")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddError("Unable to get spatial analyst extension")
    sys.exit(0)
'''
HU formula if Max Temp > Max Crop Temp HU formula if Max Temp < Max Crop Temp
1   85  (85-(MaxTemp - 85) + MinTemp)/2) - 45   ((MaxTemp + MinTemp)/2) - 45
2   75  (75-(MaxTemp - 75) + MinTemp)/2) - 40   ((MaxTemp + MinTemp)/2) - 40
3   75  (75-(MaxTemp - 75) + MinTemp)/2) - 45   ((MaxTemp + MinTemp)/2) - 45
4   80  (80-(MaxTemp - 80) + MinTemp)/2) - 50   ((MaxTemp + MinTemp)/2) - 50
5   95  (95-(MaxTemp - 95) + MinTemp)/2) - 50   ((MaxTemp + MinTemp)/2) - 50
6   90  (90-(MaxTemp - 90) + MinTemp)/2) - 50   ((MaxTemp + MinTemp)/2) - 50
7   90  (90-(MaxTemp - 90) + MinTemp)/2) - 60   ((MaxTemp + MinTemp)/2) - 60
8   80  (80-(MaxTemp - 80) + MinTemp)/2) - 65   ((MaxTemp + MinTemp)/2) - 65
9   95  (95-(MaxTemp - 95) + MinTemp)/2) - 65   ((MaxTemp + MinTemp)/2) - 65
'''
DIR = os.path.dirname(__file__)
xlsFilePath = os.path.join(DIR, '..', 'resources', 'xls')
INPUT_NORMAL_HEATUNIT_RAST_MAX_FOLDER = os.path.join(DIR, '..', 'resources',      'input', 'max')
INPUT_NORMAL_HEATUNIT_RAST_MIN_FOLDER = os.path.join(DIR, '..', 'resources', 'input', 'min')
OUTPUT = os.path.join(DIR, '..', 'resources', 'output')

for day in range(0, 367):
    maxNormName = 'Nor_Max_{0}.tif'.format(str(day).zfill(3))
    minNormName = 'Nor_Min_{0}.tif'.format(str(day).zfill(3))
    maxNormFile = os.path.join(INPUT_NORMAL_HEATUNIT_RAST_MAX_FOLDER, maxNormName)
     minNormFile = os.path.join(INPUT_NORMAL_HEATUNIT_RAST_MIN_FOLDER, minNormName)
    categories = [(maxNormFile, minNormFile, 85, 45, str(1)), (maxNormFile, minNormFile, 75, 40, str(2)),
              (maxNormFile, minNormFile, 75, 45, str(3)), (maxNormFile, minNormFile, 80, 50, str(4)),
              (maxNormFile, minNormFile, 95, 50, str(5)), (maxNormFile, minNormFile, 90, 50, str(6)),
              (maxNormFile, minNormFile, 90, 60, str(7)),
              (maxNormFile, minNormFile, 80, 65, str(8)),
              (maxNormFile, minNormFile, 95, 65, str(9))]

for category in categories:
    print "Con(\"%{0}%\" > {2}, ((({2}-(\"%{0}%\"-{2}))+\"%{1}%\")/2)-{3}, ((\"%{0}%\" + \"%{1}%\")/2)-{3})".format(
        *category)

    output = 'Nor_HU{0}_{1}.tif'.format(category[4], str(day).zfill(3))
    outputLocation = os.path.join(OUTPUT, output)

    arcpy.gp.RasterCalculator_sa(
        "Con(\"%{0}%\" > {2}, ((({2}-(\"%{0}%\"-{2}))+\"%{1}%\")/2)-{3}, ((\"%{0}%\" + \"%{1}%\")/2)-{3})".format(
            *category),
        outputLocation)

When I run this, I get a type error pointing to the last line where I am using a con statement in raster calculator. My console output is:

Con("%C:\Users\uma.bhandaram\Desktop\Peak\peakseason\normalCalc..\resources\input\max\Nor_Max_000.tif%" > 85, (((85-("%C:\Users\uma.bhandaram\Desktop\Peak\peakseason\normalCalc..\resources\input\max\Nor_Max_000.tif%"-85))+"%C:\Users\uma.bhandaram\Desktop\Peak\peakseason\normalCalc..\resources\input\min\Nor_Min_000.tif%")/2)-45, (("%C:\Users\uma.bhandaram\Desktop\Peak\peakseason\normalCalc..\resources\input\max\Nor_Max_000.tif%" + "%C:\Users\uma.bhandaram\Desktop\Peak\peakseason\normalCalc..\resources\input\min\Nor_Min_000.tif%")/2)-45)

This looks right to me but clearly it's not. Where am I going wrong?


texts:
I am trying to use Python with Grass GIS.
To get started I try to read in a raster TIFF.

The import of grass.script works fine:

import grass.script as g

But if I am trying to run a command I get the following error

sloperaster = 'slopet.tif'

a = g.run_command('r.in.gdal',
                      input = sloperaster,
                      output = 'test')

>>>'r.in.gdal' is not recognized as an internal or external command,
operable program or batch file.

Any ideas what is causing that and how to resolve it?

I am on Windows 7 (64bit), Python 2.7.5 (32bit), Grass GIS 6.4.3 (32bit)


texts: I am looking for a solution where I could group points that are nearest to each other in the given area. Having points lon/lat I would like to know how many groups with points closest to each other or by determining the number of points to find out how many points groups are located in the area. Did anyone perform such an analysis and could share the solution?


texts: How to manage unique ID's for altered polygon

texts: I'm working with administrative boundary polygons and am wondering when to keep the existing unique ID, and when to assign a new unique ID to a polygon. For example, if District A merges with District B and is still called District A, do I keep the original ID for District A, or create a new ID? What about when District A keeps the same geometry, but is renamed to District B? Or when the boundaries are slightly changed? Is there a standard way to manage unique ID's?


texts: Move shapefile from one folder to another using python

texts: I'm working with shapefiles and using the ArcPy library. My script is creating a new shapefile based on a join and I want to know how to move the newly joined-shapefile from the environment location to a new folder. Any ideas on how to do this?


texts: QGIS Plugin development version compatibility

texts:
At this time latest QGIS version is 2.4 , QT Version is 5.3.1 and PyQT5.Are these version compatible with each other ?

Once QGIS 2.4 is installed , how does one know which version of QT and PyQT(and Python too..) are required to develop Python,C++ plugins and C++ standalone applications ?


texts: I have 46 MODIS NBAR images of the same area (time-series). How can I get seven (1 for each band) multilayer (46 layer) images, for the further FastICA analysis? According to the results of the paper,  it should all result in matrix A of 46*10 in size and 10 independent source images. Input to FastICA : NDVI or SR.


texts: Georeferencing multispectral data using ArcGIS Desktop?

texts:
I often have to georeference Landsat data image to map (often raw data without geometric corrections). The way I am doing this is to open up one channel and set my GCPs and save them into a txt file. After that I open each individual channel and am loading the GCPs from the txt file to each image and I am rectifying each channel.

Is there a faster way to do this so that I don't have to load and georeference each channel separately?

I'm using ArcGIS Desktop 10.0.


texts: Align two georreferenced rasters which have different resolution in QGIS

texts:
I have two rasters in .TIFF format which have different resolution, different information and extent but georreferenced in the same coordinate system .

I would like to know if it possible to align the pixels of this two rasters to have them in the same position.

I was try with the Align raster tool in QGIS but it doesn´t works well.


texts: Getting QGIS working with Arcgis feature server - crs issue?

texts:
I am having difficulty getting the feature service below to line up properly with other layers.  It renders fine by itself.  How do I determine what crs this service is using?

I am adding this layer with the ArcGIS REST API Connector plugin.

http://services2.arcgis.com/8krRUWgifzA4cgL3/ArcGIS/rest/services/Damage1995Blowdown/FeatureServer/0


texts:
I am trying to merge  features (polygons) from 7 subregions in ArcGIS (v10.x). Merge goes fine but when I do any other process on it (e.g. Adding a field, calculations etc), the merged file is corrupted. When I want to view the shapefile I get a message "Open Failed ... Number of shapes do not match number of table records"

I tried to check geometry of the merged file (before any process was launched) and found 4 self-intersections. I tried to repair geometry to fix them but it did not help!

The version I am using is 10.2.0.3348. Unfortunately the files are heavy (each over 2GB, there are 7) to upload. Besides, they are from the Provincial Ministry and I am not sure if I could put them in the public domain.

I received the data as a geodatabase (certain forest operations in individual management units of a region). I need to extract a few variables from this data to join with another large dataset at regional scale for further analyses (décisions are at regional scale). I had earlier scripted the entire process in Python (which was working fine on several other regions). Since my script was always looking for shapefiles (and moreso I am unsure of working with geodatabases directly), I extracted individual shapefiles.

The initial process is supposed to merge the shapefiles into a region, add a few fields and calculate them. The merge goes through but no other process on it (e.g. add field).

@Polygeo, the process works on a subset (perhaps the subset did not have self-intersection issues)?

I displayed all the data in ArcMap, exported them as individual shapefiles, merged, made a new copy (through export) of the merge, and still have the same problem.


texts: Randomly extract value in time dimension at point location from NetCDF?

texts:
I have a set of point locations (point feature dataset) and a netCDF file containing 30 years of daily precipitation data.

My goal is to randomly select a precipitation value in the time dimension at each point location and write it to a table. For raster datasets, this can be accomplished by using the extract values to points tool or the extract multi values to points tool within ArcMap.

Is there a way to modify the python code for any of these tools to accomplish a random sample in the time dimension, or does anybody know of any other feasible approaches?


texts: Changing tmap custom icon/symbol size within legend

texts:
Continuing a question from here: Adding tmap_icon to legend in r
The answer provided shows how to add a custom PNG image to a tmap in R, however, I now looking for a way to increase the symbol size within the tmap legend. See below:
library(tmap)
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
library(tidyverse)

tdir=tempdir()
Stateshpurl= "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_5m.zip"
if(file.exists(paste(tdir,"/cb_2018_us_county_5m/cb_2018_us_county_5m.shp",sep=""))==FALSE){
  download.file(Stateshpurl, destfile = file.path(tdir,"State_Boundaries.zip"))
  unzip(file.path(tdir,"State_Boundaries.zip"),exdir = tdir)}
US_states=read_sf(paste(tdir,"/cb_2018_us_state_5m.shp",sep=""))

US_states = st_transform(US_states, 5070)

L48 = US_states %>%
  filter(STATEFP != "02") %>%
  filter(STATEFP != "15") %>%
  filter(STATEFP != "60") %>%
  filter(STATEFP != "64") %>%
  filter(STATEFP != "66") %>%
  filter(STATEFP != "68") %>%
  filter(STATEFP != "69") %>%
  filter(STATEFP != "70") %>%
  filter(STATEFP != "72") %>%
  filter(STATEFP != "74") %>%
  filter(STATEFP != "78")

icon_url="http://29.media.tumblr.com/tumblr_m0q2g8mhGK1qk6uvyo1_500.png"
download.file(icon_url, destfile = file.path(tdir,"Image.png"),mode = 'wb')

L48$icon = "my label"
Icon=tmap_icons(file=paste(tdir,"/Image.png",sep=""))

tm_shape(L48)+
  tm_borders()+
  tm_symbols(shape = "icon", shapes = Icon,
             title.shape = "My icon")

Notably, the size =  argument only changes the size of the icon within the map space itself.


texts:
Given a set of 2D points ( or vertices if you want ), the values stored for each point are:

coordinate on x
coordinate on y
an enum or a value from a list in general

for the sake of this example let's say that this last value is the name of a fruit, so the record in the data structure for each point looks like

"34.3" "845.2" "apple"

The set of fruits considered is finite and well defined at the moment the algorithm comes in to play, with this I mean that I can easily list all the possible variations for that 3rd value, if I need to.

Since in my model that 3rd value is assigned per-face, I don't really care if my algorithm acts on faces or on vertices; I should also note that the number of vertices for each face are not equals, so with the term "face" I mean a generic geometry that can be a triangle, a quad, and so on, but no matter how many vertices a face owns, the property is given per-face.

Now the real problem: I need to do some clustering based on that 3rd arbitrary value.

I don't have any pattern here, I have no guarantee on anything, for what I know it's possible to have the entire map composed of single faces with a different 3rd value, they can be convex or concave, the only thing is, I have this kind of records in input, and as output I need a set of areas grouped according to this 3rd value that they have in common.
Offcourse if 2 faces share the same 3rd value but they are not adjacent ( don't share at least 1 edge/verts ) they should be considered as 2 separated areas.

What is an algorithm for this kind of clustering ?

Since the dataset is relatively small the time-complexity of the algorithm doesn't really matters that much.

PS
If needed I can simplify every face down to a triangle; I am also not sure if this problem can be solved with a Convex Hull algorithm.


texts:
I am trying to plot district boundaries on a map of Nepal.  For that I  downloaded a shapefile of districts in Nepal from this site:

  http://www.arcgis.com/home/item.html?id=d4ad6faa446b42bd9f1d42407db7b5ba

The folder name was  NP_75DWGS84 and all the files within the folder had the files containing the same name. Then I used the following code:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap()
map.drawmapboundary()
map.drawcoastlines()
map.readshapefile('/home/mala/NP_75DWGS84/NP_75DWGS84','NP_75DWGS84',linewidth=1.0,drawbounds = False)
plt.show()

It gives me total plot of world, but, it is not plotting shapefiles in the graph. There is also no error. What is wrong here?


texts: Automatically adding dimensions to polygons in QGIS

texts:
Is it possible to show :

lot #
area
dimensions of the polygon's lines
highlighted polygon

with the only input being the lot #. The result should look like this:


texts:
I have some updated waypoint in one layer that correspond to waypoints in another (the 'Master') and I have a script that is supposed to update the geometry in the 'master' from the update layer.

layers = QgsMapLayerRegistry.instance().mapLayers()

# find the update and master layers

master = update = None
for name, layer in layers.iteritems():
    if 'wp_master' == layer.name():
        master = layer
    elif 'changed'  == layer.name():
        update = layer

update_features = {}  # indexed by feature attribute name
for f in update.getFeatures():
    update_features[f.attribute('name')] = f

# now search the master layer for the ones that have changed

master.startEditing()
for m in master.getFeatures():
    n = m.attribute('name')
    if update_features.get(n, 'none') != 'none':
        u = update_features[n]
        old = m.geometry().asPoint()
        new = u.geometry().asPoint()
        print 'master', n, old.x(), old.y()
        print 'new', n, new.x(), new.y()
        # if not m.geometry().equals( u.geometry()) :
        print 'updating', n
        print 'time', time.strftime("%Y-%m-%d %H:%M:%S")
        m['updated'] = time.strftime("%Y-%m-%d %H:%M:%S")
        master.dataProvider().changeGeometryValues({m.id(): u.geometry()})
        master.updateFeature(m)
        master.commitChanges()
        new = m.geometry().asPoint()
        print 'nex x,y ', new.x(), new.y()

gives this output:

Running mod-test.py in: /Users/rful011/src/qgis/scripts/bin
master T J B03 3 1768944.68311 5947703.25154
new T J B03 3 1768943.24096 5947702.72496
updating T J B03 3
time 2018-02-25 11:00:38
new x,y  1768944.68311 5947703.25154
Completed script: mod-test.py

We note that the printed geometry of the point in master is unchange

What am I missing and is there a better way of doing this without iterating over the master layer?


texts: Integrating SAGA 7.8.1 into QGIS 3.16 (OS Windows 10)

texts:
I want to calculate the PV potential of roofing using QGIS, similar to this post:
How to calculate total PV potential of a roof using QGIS?.
However, my SAGA version in QGIS does not show me the Potential Incoming Solar Radiation algorithm.
I have downloaded a standalone SAGA version 7.8.1.
Can I integrate the latest/ update the current SAGA version that is used by my QGIS?
The OS I am using is Windows 10.
I am completely new to SAGA in QGIS and have no previous experience in using it.


texts: Too many features selected with the ArcMap Select by Polygon tool

texts:
I have ArcMap 10.1 SP1. My features are stored in Oracle 11g database in SDO_GEOMETRY format.

I have imported a layer of accidents location in ArcMap. So when I try to select all features included in a polygon (using the Select by Polygon), then too many features are selected. I mean some feature outside my polygon are part of selected feature.

I draw a triangle on my map and I expect that all features included into the triangle will be selected. The result I obtain is all features included in the rectangle envelope of the triangle that are selected, which is obvious wrong.

Is someone can help me on this one?

I was in SDE_GEOMETRY before and I did not have the problem when I tried to select by polygon.


texts: Import Esri linear referencing system into Oracle Spatial

texts: I am wondering if it's possible to import a linear referencing system (routes with measures) into Oracle Spatial?


texts: Implementing a category theme using SharpMap.Rendering.Thematics.CategoryTheme?

texts: GitHub page for SharpMap has a link for related solution for CategoryTheme without any sample code.  Can someone provide sample code using the SharpMap.Rendering.Thematics.CategoryTheme for a shapefile using SharpMap?


texts: Coding Select Polygon and Make New Layer

texts:
I have a large grid file, and I am looking to extract each individual polygon that makes up the grid into a new layer. Essentially, I would like to replicate selecting a layer by mouse click, then "Create Layer from Selected Features".

I am familiar with coding in ArcPy and the Python environment as a whole, I suspect the best approach is to script this process there. My grid table does have each cell attached to a unique number, should this prove useful.

I've written the following code to attempt this:

import arcpy
import sys
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/Users/Eric/Documents/"

#Set Local Variables
in_grid = "fishnet_clipped.shp"
out_feature = "polygrid2.shp"
where_clause = ' "FID" = \'2\' '

# Executing Select
arcpy.Select_analysis(in_grid, out_feature, where_clause)


texts: Leaflet - map.fitBounds() / map.getBoundsZoom() BROKEN after resizing map?

texts:
https://jsfiddle.net/MaxHeadroom/s6h40xrn/

If you resize the map size, the zoom level should be adapted to always show the orange rectangle at maximum possible size.

But fitBounds seems to be broken after resizing the map.
Same problem with map.getBoundsZoom().

To reproduce the buggy behaviour, start with a small map size - increase it and decrease it then.
Check the console output too.

Or is there some a mistake in my source code?


texts:
I would like to use the QGIS API in my own program to display maps ....

So I thought I would find some tutorial here :
http://hub.qgis.org/wiki/quantum-gis/Writing_custom_Applications_with_QGIS_API

But the links are directed towards the qgis blog (http://blog.qgis.org/?q=node/10), which seems to be "down" ....

Where can I find such tutorials ?

Best regards.


texts: Getting values from another column, from different row but from same table using QGIS?

texts:
I am using QGIS 3.4.

I have valves on a water pipeline that are in either 'on' or 'off' ('tilstand' - 'on' / 'off') condition.
In another column, there is the valve, that is physically placed before the valve in the column "no".
I manually fill in column "v_hl_ventil-1" but want column "tilstand-1" to automatically retrieve data from column "tilstand" for the corresponding “nr” and change automatically, depending on condition in “tilstand”.


texts:
How can I edit a QGIS project file so that it would open, save and edit data from given folder and save the results for different directory? Kind of like MapInfo's Wor. file.

I read from QGIS manual that the project file is saved in a XML format, but for some reason I can't open it.


texts: Add feature layer to display with in a script tool for desktop mapping 'App'

texts:
The point of this question is to find a way to programatically change the map extent based on user input.  My method involves creating a Feature Layer & Selection, but if you have another (different) solution that gives the desired result, I'd love to hear it!

Now that I'm comfortable with arcpy and I've become aquainted with the mapping module I saw a niche.  I can write start-to-finish desktop mapping applications to accomplish a specific task.  This enables anyone in our office to run the 'app' without needing several years of using ArcMap to accomplish a repeatable task for a different area of interest (ie. zip-code, tax-parcel, county, etc.)

I've written a few of these scripts and associated each of them with their own MXD to give the user a GUI.  When they input some sort of attribute value (watershed number, land-parcel number, county, etc.) the script makes a feature layer of that polygon and uses the zoom to selected feature method to change the extent.

The problem I'm having once I've finished debugging the 'app' through ArcMap's Python Window I create a script tool that gives the user an input box to enter in the required arguments for the program to run.

For whatever reason, ArcMap defaults to NOT adding the feature layer to the map when it reaches the line of code where it is created.  Therefore, my script cannot use the zoom to selected feature method to change the extent.  The rest of the script runs along just fine, but the zoom extent is still at the state level as opposed to a sub-watershed level.

Code excerpt:

try: # Clip wetlands layer to extent of user dictated parcel, calculate staticstics
    newWetland = outWetlands + "HUC_" + userHucID # variable that stores clipped wetland output location & name
    arcpy.MakeFeatureLayer_management(sheds, "currentShed", clipClause) # Create feature layer of HUC (from user input)
    arcpy.SelectLayerByAttribute_management("currentShed", "NEW_SELECTION", clipClause) # Select current HUC
    df.zoomToSelectedFeatures() # Zoom to extent of selected HUC
    arcpy.RefreshActiveView()   # Refresh view of mxd before exporting to PNG

    arcpy.Clip_analysis(wetlands, "currentShed", newWetland) # Clip dissolved wetlands layer to current HUC
    arcpy.ApplySymbologyFromLayer_management(newWetland, layers + "WetlandSymbology.lyr") # Apply wetlands symbology
    print "Finished creating wetland file for HUC #: " + userHucID + "."
except:
    print arcpy.GetMessages()
finally: # delete the feature layer if it exists if an error is present or not
    if arcpy.Exists("currentShed"):
        arcpy.Delete_management("currentShed")

Script run in the Python Window of ArcMap:

SAME script run as Script-tool in custom toolbox:

Script tool Properties window:
This was suggested to me by an ESRI rep -- add a feature layer as a derived, output parameter, then use the SetParameter(index, featureLayer) to get the desired layer added to display.

The parameter is derived (it isn't given by the user, but is derived by one of the user inputs), and set to output in an attempt to add it to the map.

Attempted Fixes:

ArcMap Toolbar -> Geoprocessing -> Geoprocessing Options -> 'Add results of geoprocessing operations to the display'
Add Feature Layer output parameter in script tool properties paired with arcpy.SetParameterAsText method
arcpy.mapping.AddLayer(df, layerPath)
df.extent = layer.getExtent()


texts: Modelling fiber optic network using free solution in GIS?

texts:
I am a beginner to GIS software and currently I have been given a task to model a telecom network. Most of the solutions that I came across are proprietary.

I want to know which open source software can I use to show the connectivity between multiple fibers (FTTX network)?

I do not need a software to plan and design a network, I just have to make a replica of the network that I already have in KMLs.


texts:
So, I have a few rasters I need to use for processing, but the task requires them to be in a different SRS than they currently are. GDAL Warp would be the go-to solution, however, it saves the transformed raster in a new file, and for project reasons I cannot have it save a transformed file anywhere, it has to be done in memory. I thought of two ways to do it, but encountered problems in both.

  Alternative #1: use the MEM driver.

  This would be ideal, but GDALDriver::CreateCopy() however does not
  have a reprojection parameter, and I'm unable to create an empty
  memory raster (and pass that to Warp) because GDALDriver::Create()
  expects input of the number of rows and cols, and I have no idea how
  many there'll be in the transformed raster before I transform it.

or

  Alternative #2: build a VRT.

  I can very easily get the desired data from a transformed vrt, and
  that would actually be my preferred solution if I was able to save it
  anywhere. But I can't. But then I thought, a vrt is just an xml, might
  there be a way to just generate the vrt in memory, instead of saving
  it to disk beforehand?

I'm not really tieing this problem to one alternative or the other, just showing what I tried thus far. A working solution for either would be fine, as would be an entirely different solution I haven't thought of.


texts: Why does GDAL Polygonize flip features vertically?

texts:
I am trying to convert a 1 band jpg raster into polygons.
It is converting the data and I get a polygon shape file. However, the output data is vertically flipped.
Input data have no reference system.

I am using the code in c++.


texts:
Starting with an .img file how to I convert a specific (x,y) coordinate into (lon,lat) with gdal. I don't know what my input format is, but if I type gdalinfo topo.img it tells me the corner coordinates of the map, so it ought to be possible to get the coordinates for any of the grid points without further input. gdallocationinfo topo.img [x] [y] gives the elevation only.

The beginning of the gdalinfo output for the input file looks like this:

Driver: HFA/Erdas Imagine Images (.img)
Files: bigislandDTM.img
       bigislandDTM.ige
       bigislandDTM.rrd
Size is 32609, 35573
Coordinate System is:
PROJCS["NAD_1983_UTM_Zone_5N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-153],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["EPSG","26905"]]
Origin = (170717.059624216170050,2247879.339183096773922)
Pixel Size = (4.497066284309543,-4.497066284309542)
Corner Coordinates:
Upper Left  (  170717.060, 2247879.339) (156d 9' 9.63"W, 20d18' 1.76"N)
Lower Left  (  170717.060, 2087905.200) (156d 7'28.87"W, 18d51'25.08"N)
Upper Right (  317361.894, 2247879.339) (154d44'57.59"W, 20d19'12.58"N)
Lower Right (  317361.894, 2087905.200) (154d44' 1.61"W, 18d52'30.47"N)
Center      (  244039.477, 2167892.270) (155d26'24.29"W, 19d35'22.45"N)


texts: Python OGR Geometry changes Floating Points

texts:
I'm creating a Polygon in OGR and adding coordinates to it in Python. It is a must that my coordinates aren't allowed to change when adding, however they do and I don't understand why. Here is my example:

import ogr
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(12.5433, 54.10064)
ring.AddPoint(12.5433026066, 54.1006582464)
ring.AddPoint(12.543312311, 54.1006489082)
ring.AddPoint(12.5433093954, 54.1006320625)
ring.AddPoint(12.5433, 54.10064)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

polygon = poly.ExportToWkt()

Now if I print the polygon, I receive:

'POLYGON ((12.5433 54.10064 0,
           12.543302606599999 54.100658246400002 0,
           12.543312311 54.1006489082 0,
           12.5433093954 54.100632062499997 0,
           12.5433 54.10064 0))'

Clearly floating points have been added to some coordinates and I can't understand why.

Example: (12.5433026066, 54.1006582464) to (12.543302606599999 54.100658246400002)


texts:
I'm creating a territory map using QGIS that is based on both county and ZIP code. I need a shapefile for all county/ZIP combinations in Arizona (would look like a ZIP code shapefile except that all ZIP codes that are in more than one county would be split by county lines into multiple polygons).

Is there an easy way to do this?


texts: Drawing shapes with highest possible decimal accuracy

texts:
I am trying to read input of lat and long to draw the shape based on condition.The input is standard array of lat and long with 15 decimal accuracy. I am realizing - I am loosing some decimal accuracy (the point was somewhat offsite and complex polygons were never drawn).
The samples are here:

point: [[[-81.00392112682097,46.50255058560673]]]

the plotted point is: -81.003906    46.502686 and I am generating PDF maps at 1:3000 scale so it is a big deal.

polygon: [[[-79.37485174735863,43.6430783790403],[-79.37361793121188,43.64341222822559],[-79.37323169311376,43.64261253986735],[-79.374186559523,43.64234079868168],[-79.37412218650664,43.64220880966229],[-79.37454061111293,43.64210787668709],[-79.37471227248987,43.64248831547784],[-79.37463717063745,43.64251160757044],[-79.37481956085045,43.64285322389105],[-79.37470154365381,43.642884279823875],[-79.37485174735863,43.6430783790403]]]

I used eval to read input as list, but probably that was not a good idea. Not sure how and where should I cast the input as decimal?

OrderCoord =  eval(arcpy.GetParameterAsText(0))
point = arcpy.Point()
array = arcpy.Array()
 # A list that will hold each of the Polygon objects
featureList = []
for feature in OrderCoord:
    # For each coordinate pair, set the x,y properties and add to the
    #  Array object.
    #
    for coordPair in feature:
        point.X = coordPair[0]
        point.Y = coordPair[1]
        array.add(point)

    if OrderType.lower()== 'point':
    feaERIS = arcpy.Multipoint(array)
    elif OrderType.lower() =='polyline':
    feaERIS  = arcpy.Polyline(array)
    else :
    feaERIS = arcpy.Polygon(array)
    # Clear the array for future use
    array.removeAll()

    # Append to the list of Polygon objects
    featureList.append(feaERIS)

I added the screenshot of these points at 2000 scales. We have a Flex application where customer can draw and order a geometry. All the data sits in ORACLE and I am working on a GP Tool to print final PDF.


texts: esriPFS(Picture Filled Symbol) imageData is not displayed on Leaflet Map

texts:
import * as L from "leaflet";
import * as esri from  "esri-leaflet"
import "esri-leaflet-renderers";

let _layer = esri.featureLayer({
url: 'https://services.arcgis.com/aA3snZwJfFkVyDuP/ArcGIS/rest/services/airport_districts_and_spheres/FeatureServer/1',
simplifyFactor: 1
});

this.map.addLayer(_layer);

Since Esri-leaflet-renderer does not support "symbol-type: esriPFS", is there any way we can render esriPFS (Picture Filled Symbol) on Leaflet map?


texts:
My script is intersecting lines with polygons. It's a long process since there are more than 3000 lines and more than 500000 polygons. I executed from PyScripter:

# Import
import arcpy
import time

# Set envvironment
arcpy.env.workspace = r"E:\DensityMaps\DensityMapsTest1.gdb"
arcpy.env.overwriteOutput = True

# Set timer
from datetime import datetime
startTime = datetime.now()

# Set local variables
inFeatures = [r"E:\DensityMaps\DensityMapsTest.gdb\Grid1km_Clip", "JanuaryLines2"]
outFeatures = "JanuaryLinesIntersect"
outType = "LINE"

# Make lines
arcpy.Intersect_analysis(inFeatures, outFeatures, "", "", outType)

#Print end time
print "Finished "+str(datetime.now() - startTime)

My question is: is there a way to make the CPU work at 100%? It's running at 25% all the time. I guess that the script would run faster if the processor was at 100%. Wrong guess?

My machine is:

  Windows Server 2012 R2 Standard
  Processor: Intel Xeon CPU E5-2630 0 @2.30 GHz 2.29 GHz
  Installed memory: 31,6 GB
  System type: 64-bit Operating System, x64-based processor


texts:
There does not seem to be a command to get the number of layers in the map, or am I missing it?

One way to count the layers would be:

Object.keys(map._layers).length

Or is there another preferred way? For instance, I would be interested in finding out how many vectors or how many basemaps have been added to the map, and I am not sure how to achieve that via the above snippet.


texts: Is panchromatic usable for supervised classification schemes in R or QGIS?

texts:
I am using Worldview-2 and Worldview-3 panchromatic images from 2019. I am looking at discrete topographical features which are easily delineated as polygon ROIs. However, I only have this panchromatic imagery to work with.

The goal for this project is to delineate tindividual, discrete topographic features as polygons in QGIS, use them as training files for classification, and use a classification method to identify topography within the image in R or in QGIS.

I have tried classifying the images using tools such as SCP, and functions in the Orfeo Toolbox. Although, it seems all of these tools require multiscpectral imagery.

Is image classification only possible with multispectral imagery, or can images be classified using panchromatic imagery? What tools allow for the latter?


texts: Searching and downloading OpenStreetMap data in QGIS

texts:
My OpenLayers plugin is under web and not under the plugin and I can't find add OpenStreetMap  layer. I went ahead and opened web > open layers plugin > openstreetmap but this is the error I keep getting.

  2018-04-10T13:40:26   1   Traceback (most recent call last): File
  "C:/Users/RUKEVWE/.qgis2/python/plugins\openlayers_plugin\weblayers\weblayer.py",
  line 100, in addLayer self._addLayerCallback(self) File
  "C:/Users/RUKEVWE/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py",
  line 182, in addLayer self.setMapCrs(coordRefSys) File
  "C:/Users/RUKEVWE/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py",
  line 231, in setMapCrs extMap = coordTrans.transform(extMap,
  QgsCoordinateTransform.ForwardTransform) QgsCsException: forward
  transform of (-8.208501, -2.390302) PROJ.4: +proj=longlat +datum=WGS84
  +no_defs +ellps=WGS84 +towgs84=0,0,0 +to +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 Error: latitude or longitude exceeded limits


texts: How to use GDAL translate tool in python script QGIS

texts:
I'm trying to use gdal_translate tool in a python script I'm writing in QGIS 2.18.3. I use the translate tool because I want to change all the zero values (in all 8 bands of my raster) to nodata. I'm able to do so manually by using the tool from the toolbox. When I run my code in QGIS script editor the processing window disappear after a second, I don't get any error message but no new file is created. I think it is something to do with how this tool understands what is raster, but I've successfully used other gdal tools (merge, raster calculator) in python script using the same method, so I don't understand why it does not work in the translate tool.

This is a simple code I'm trying to use:

##last_try=name
##Select_directory=Folder

import glob, os
from PyQt4.QtCore import QFileInfo
from qgis.core import *

for raster in glob.glob("*.tif"):
    fileInfo = QFileInfo(raster)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(raster, baseName)
    trans = os.path.join(Select_directory,"image_suffix.tif")
    outputs_GDALOGRTRANSLATE_1=processing.runalg('gdalogr:translate', rlayer,100.0,True,'0',0,'',None,False,5,4,75.0,6.0,1.0,False,0,False,None,trans)

As I mentioned, this type of code works great with other gdal tools.

What's wrong with my code?


texts: I want to query features under [Layer Properties] -> [Source] -> [Query Builder] but the query builder button is greyed out. The layer I'm using is from a shapefile with multiple demographic fields by US county. I'm using QGIS 3.16.


texts: How do i get ArcObjects not to ignore a layer definition that has been set?

texts:
I am creating a simple text search application, where the user can enter an address search and the result appear in a listbox, which they can then select it and zoom to the address.

Now the problem I am having is that I have set  layer definition, within the layer I am working with to only show Residential addresses ( and this works, as i can only see residential addresses in the attribute table), however when I enter a street (in the app), that contains not only residential address, the records appear in my list box when though I have set the layer definition.

Is this behavior normal ? and what is the best way to get round this?

Thanks


texts: How do I manipulate a line layer to have straight features instead of curves?

texts:
I am working with a layer (in yellow) that i want to make follow the red line.  I have tried:

Snapping - this results in geometry issues where a line intersects itself
Smooth - has no effect

Addtionally the red layer does not cover the same extent as the yellow layer and so the snap only moved parts of the yellow layer.

Ultimately the yellow layer needs to follow the points in the same way the red layer does.

I am working out of QGIS and only have access to QGIS (not ArcMap etc)


texts:
QGIS V-3.4.5, Python 3.7.0

from qgis.gui import *
from qgis.core import *

import sip
sip.setapi('QString', 2)

from PyQt5.QtCore import * # QSettings, QTranslator, qVersion, QCoreApplication

from PyQt5.QtGui import * #QIcon
from PyQt5.QtWidgets import * #QAction

from PyQt5.QtXml import QDomDocument

import os.path
import os

...

 actionSelect = QAction( QString("Select"), self)

             NameError: name 'QString' is not defined

NameError: name 'QString' is not defined

WHY???


texts: Drawing polygons that share boundaries using Google My Maps

texts:
I am relatively new to GIS and I am facing a problem to maintain a layer of polygons that share boundaries. I use google my maps to edit the polygons but they are full of gaps and overlaps.
How can I create polygons that share boundaries and have no gaps and overlaps?


texts: How to color markers based on polygon area?

texts:
I have a map with around 317 markers and each marker is the center point of a polygon. The polygonal area varies throughout, so I want to represent the markers in such a way that the polygon with the least area has the lightest shade of the marker color and the marker representing the largest area has the strongest shade of the same marker color.

Can I dynamically produce such colored icons on google maps api??


texts: how does arcmap calculate scalebar inside a WGS84 layout?

texts:
ArcMap lets you put a scalebar inside a layout referncing a WGS84 data frame.
How is the scale calculated? Is it a "mean" scale? And, in case it's right, what "mean" scale could mean?


texts: Polygon construction Do's and don'ts Advice?

texts:
I am working on two polygons. One is Forest layer that has the ownership and I have another polygon that has the fire perimeter. I need to build polygon construction to cut out some of them like we need to take out small area like non-forest.

For example, I have the fire perimeter and I need to fire polygon overlay the forest ownership but there are some area that doesn't overlay right. The Line in between non-forest and Forest ownership and I don't want the fire perimeter to over the non-forest but just the forest ownership. Some overlay the private areas that not need to be shown.

I do not know if this is the right method to clean it up after I contsturct the polygons.

Help ?

You can read this from ArcGIS 10.2 help about it

http://resources.arcgis.com/en/help/main/10.2/index.html#/Constructing_polygons_from_the_shapes_of_other_features/01m700000012000000/


texts: Blank Excel fields creating problem in Creating TableView using ArcPy?

texts:
I was about to take input from an excel and using conversion tool in arcpy I have to convert the excel to table in a gdb. But when I used this Excel to table tool, I found that there are columns named 'F34', 'F35' etc which in excel does not hold anything.

How do I either remove those columns from the table or avoid those columns while doing conversions?


texts: How to calculate polygon areas and perimeters using QGIS?

texts:
Using QGIS, if I load a polygon layer and make a new column, how do I fill that column with the:

area of each polygon and another column with
length of the perimeter of each polygon?


texts: How do I set the Coordinate Reference System for a NetCDF file through a configuration file in QGIS?

texts:
I have downloaded some NetCDF files (.nc) from the following link: http://fizz.phys.dal.ca/~atmos/martin/?page_id=140. However, when I open these files in QGIS, they do not align with my vector layer.

To realign them, I currently convert the NetCDF filed to GeoTIFF files and then specify the following CRS through a TFW file:

0.0100000000
 0.0000000000
 0.0000000000
-0.0100000000
-179.9950000000
 69.9950000000

Is there a way for me to set the CRS for NetCDF files directly through  configuration files, without converting them to GeoTIFF?


texts: Elevation profile sample point algorithm

texts: In generating an elevation profile, samples from the underlying raster must be taken. If the objective is to create a profile of maximum quality, how should the locations of sample points be determined? What approach should be taken to interpolating those sample points?


texts: Inserting Geometry into PostGIS Table Violates Check Constraint "enforce_geotype_geom"

texts:
I have imported a shapefile containing point layer into postgis, and the SRID i have chosen was 4326. Now while inserting the values into the table I am getting a problem. My intension is to create a circle around that point using ST_Buffer.

My Table structure is

CREATE TABLE "ankCircle"
(
  gid serial NOT NULL,
  id numeric(10,0),
  geom geometry,
  CONSTRAINT "ankCircle_pkey" PRIMARY KEY (gid),
  CONSTRAINT enforce_dims_geom CHECK (st_ndims(geom) = 2),
  CONSTRAINT enforce_geotype_geom CHECK (geometrytype(geom) = 'POINT'::text OR geom IS NULL),
  CONSTRAINT enforce_srid_geom CHECK (st_srid(geom) = 4326)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE "ankCircle" OWNER TO postgres;

Query I am using is

INSERT INTO "ankCircle"(gid, id, geom)
VALUES (1, 1, ST_Buffer(ST_GeomFromText('POINT(100 90)',4326),50, 'quad_segs=8'));

error got

ERROR:  new row for relation "ankCircle" violates check constraint "enforce_geotype_geom"

********** Error **********

ERROR: new row for relation "ankCircle" violates check constraint "enforce_geotype_geom"
SQL state: 23514


texts: Why ArcGIS JavaScript text label does not show on feature?

texts:
My working example here does not show the labeling information on features.
    var featureLayer;

    map = new Map("map", {
      basemap: "satellite",
      center: [-46.807, 32.553],
      zoom: 3
    });

    //create a feature collection for the flickr photos
    var featureCollection = {
      "layerDefinition": null,
      "featureSet": {
        "features": [
            new Graphic({
                geometry: new Point(-47.476661, 42.702002, map.spatialReference),
                attributes: {ObjectID:1 ,title: "my title", description: "my description"}
            }),
            new Graphic({
                geometry: new Point(-46.476661, 31.702002, map.spatialReference),
                attributes: {ObjectID:2 ,title: "my title", description: "my description"}
            })
        ],
        "geometryType": "esriGeometryPoint"
      }
    };
    featureCollection.layerDefinition = {
      "geometryType": "esriGeometryPoint",
      "objectIdField": "ObjectID",
      "drawingInfo": {
        "renderer": {
          "type": "simple",
          "symbol": {
            "type": "esriSMS",
             "style": "esriSMSCircle",
             "color": [216, 15, 20, 255],
             "size": 18,
            }
        }
      },
      "fields": [{
        "name": "ObjectID",
        "alias": "ObjectID",
        "type": "esriFieldTypeOID"
      }, {
        "name": "description",
        "alias": "Description",
        "type": "esriFieldTypeString"
      }, {
        "name": "title",
        "alias": "Title",
        "type": "esriFieldTypeString"
      }]
    };

    featureLayer = new FeatureLayer(featureCollection, {
      id: 'my-layer'
    });

    var labelClass = new LabelClass({
        "labelExpressionInfo": {"value": "{title} title"},
        "labelPlacement":"above-right",
        "symbol": new TextSymbol()
    });

    featureLayer.setLabelingInfo([ labelClass ])

    map.addLayers([featureLayer]);


texts:
I doing a project on Shapefile. I am trying to convert FDO Geometry to GeoAPI Geometry.

suggest me.


texts:
I'm using this vendor option to get the labels to be aligned with the
polygon :
<VendorOption name="polygonAlign">mbr</VendorOption>

As stated in the description, "If the label does not fit horizontally the
minimum bounding rectangle will be computed and a label aligned to it will
be tried out as well"
I'd like my labels to always be aligned even if it fits horizontally.
Is there a way to force this ?
I've looked at the source code but can't find the place where the Vendor Options are enforced.
Thanks


texts: How to set the viewport the bounds what I need

texts:
I need to get the world map on the screen limited to the bounds like such latlon coordinates as below. Is it possible to convert these geocoordinates to pixels and then to set these pixels to the viewport or to do this in some other way? I try to manipulate with the Rectangle2D parameters, but it is not the solution of this problem. Maybe you can give some rescue suggestions, please?

The end result of this application should be the part of the world map in the stereographic projection limited to these coordinates and saved as a picture of at least 1600 x 1200 pixels in size and at the same time with acceptable image quality. The excess territory is needed to be removed. I saved the too large image of the world map (see picture), and then I used the MarvinImage library to remove the extra territory, but at the same time the quality of the final image is greatly reduced. Therefore, I try to enlarge the the desired area of the world map, which is displayed in the JFrame, so it could be saved as a picture without additional cutting. In other words? I need to enlarge so the worl map to get the geocoordinate xy_1 in the left top corner of the JFrame, the xy_2  - in the right top corner and so on, and the size of JFrame must be close to the real size of the monitor screen.


texts: Clear markers from OSM map with openlayers

texts:
I have a bunch of Point data that refreshes every few minutes. To render these points, I have created the markers as OpenLayers.Feature.Vector Objects and added them to my Vector layer. The layer gets added to the map object, and I can see the markers when the page loads.

When the data refreshes on my local server, some points from the previous dataset may have disappeared and some new points may have been added. I am just given a large list of points to render every refresh.

Because of this, I need to clear the vector layer of all current markers and then re-render the map with the new dataset.

Also, is there a way to delete specific markers given just a lat,lng? Or is there some way to reference them and delete it from the vectory layer?

How do I clear the vector layer of all markers and how do I delete specific points?


texts: Writing custom QGIS function that works without a feature given as input

texts:
Lets say I have created a custom QGIS function, like for example
from qgis.core import *
from qgis.gui import *
import urllib.request
import urllib
import json

@qgsfunction(args='auto', group='Custom')
def jsonfromurl(url, feature, parent):
    header = {"accept":"application/json"}
    request = urllib.request.Request(url)
    response = urllib.request.urlopen(request)
    response_data = response.read()
    encoding = response.info().get_content_charset('utf-8')
    data = json.loads(response_data.decode(encoding))
    #txt = json.dumps(data) #uncomment this line to receive data as string. also change next line to: return txt
    return data

This function only works when evaluated against a feature. Otherwise it returns the preview No feature was found on this layer to evaluate the expression. and just nothing when applying it.
How can I change that function, so it also works if no feature is given. For example (used as expression) in a print layout text box or an empty layer?
Like here: Adding a textbox to a QGIS layout where I want to display the result from the request https://nominatim.openstreetmap.org/reverse?format=json&lat=51.6895477895927&lon=-0.43032246110111033&zoom=18&addressdetails=1.

Simply removing feature and/or parent seems not to work and as stated on https://howtoinqgis.wordpress.com/2017/05/20/how-to-create-custom-functions-in-qgis-using-the-function-editor/ these arguments are always needed.


texts:
I am creating markers using this method :
var markers = new OpenLayers.Layer.Markers( "Markers" );
map.addLayer(markers);
var size = new OpenLayers.Size(21,25);
var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png',size,offset);
`enter code here`markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(0,0),icon));

Now , I want to rotate the marker 90 degrees to the left (or right).
How can I do it?
I was able to find example to rotate a feature , but I am not using a feature.
Instead, I am using a marker.


texts: PS/SBAS Processing of S1 in Envi Sarscape 5.2

texts:
I'm processing a stack of s1a images in PS. I am now running the inversion first step.

In the progress bar, it is stuck at the 'preliminary candidates selection' for quite some time. I am not sure if the process has crashed or is it just taking a long time to process.

Anyone knows if this is normal?


texts: Batch/Bulk Insert ADO.Net collection in MapXtreme Table

texts:
I am working on MapXtreme 7.1 and using MICommand to insert data in a table. I have more then 10,000 records and i am inserting these using a loop. That is

Geometry = new MapInfo.Geometry.Point(coordSys, new MapInfo.Geometry.DPoint(Longitude, Latitude));
cmd.Parameters[0].Value = Geometry;
cmd.Parameters[1].Value = new MapInfo.Styles.SimpleInterior(36, System.Drawing.Color.RosyBrown);
cmd.CommandText = "Insert Into MyTable(MI_Geometry,MI_Style,latitude,longitude) values(geometry, style, latitude, longitude);
cmd.Prepare();
cmd.ExecuteNonQuery();

Here MI_Geometry and MI_Style are of types MIDbType.FeatureGeometry and MIDbType.Style respectively. The above code is working in a loop and it is inserting all of my data into that table.

The problem is, it is taking too much time to populate. Is there any way to bulk/batch insert or direct populate this table from my .Net collection through some sort of binding etc. My first requirement is to draw more then 10,000 points on Map from my collection.


texts: Redraw info window after centering map on click (ESRI Javascript API)

texts:
When I click on one of the features on my map the map gets re-centered using the following event handler:

var pointClick = projects.on("click", function(e){
    // Center map on selected feature
    map.centerAt(e.mapPoint);

    // A bunch of other stuff
});

Unfortunately the info window describing that point does not move to the new location unless the map gets redrawn (by moving to another zoom level, for example).

Is there a method to force the info window to refresh?  The information in the info window reflects the point that was most recently clicked on, but the location of the info window does not change.

Screencast:
http://www.screencast.com/t/fo1Nec0CBzk6

I can provide the code if necessary.


texts: How to properly calculate annual irradiance using QGIS

texts:
I want to calculate annual irradiance for a couple of buildings to establish how much PV energy can be potentially harvested.

Which algorithm You recommend for this purpose: "r.sun" from GRASS toolbox, or SAGA "Potential Incoming Solar Radiation" ?

Can You recommend any tutorial that explains step by step what inputs are required (i.e. not only DEM/DSM but also any atmospheric factors), what settings (or formulas) should be applied to obtain accurate results?

I am aware that this subject (or its elements) might be already discussed, but since I am a newbie, I am looking for comprehensive, holistic explanation what approach I should follow.


texts: MinimizeRestore Maximize boxes missing in Windows 10 task bar

texts:
I am running QGIS 2.18.15 on a Windows 10 OS 64-bit PC and the click symbols to close (x), minimize, restore, maximize are missing. I tried uninstalling and reloading the latest version and had same problem. I can use the keystroke commands but it is impossible to resize the window that way. In order to move the QGIS window between screens I must also use the keystroke commands - cannot drag between screens.

Everything else seems to work fine.

Am just wondering if this is a common problem and if there is a fix.


texts: Clipping geometries in PostGIS table if they are inside polygon?

texts:
I have polygons defined in KML files just like this one, which I can easily parse. Example:

<Polygon><outerBoundaryIs><LinearRing><coordinates>
-91.463940,40.372660,0
-91.483790,40.379440,0
-91.493780,40.392410,0
</coordinates></LinearRing></outerBoundaryIs></Polygon>

I would like to run queries on a postgis table and clip all lines that intersect this polygon. Segments of the lines within the polygon should be deleted. Something like this:

    ____
___|poly|___line
___|    |___line
   |____|

Here's the table structure:

    Column    |     Type
--------------+---------------
 ogc_fid      | integer
 wkb_geometry | geometry
 id           | numeric(8,0)
 ele          | numeric(12,3)

Is this possible?

What would a query look like?

I assume I need to express the polygon in WKT format and use the documented ST_Intersects method, but I'm not sure how to put this together and modify the geometries.

Can this be done with PyQGIS?

I've loaded my kml file as a QgsVectorLayer, I can access the polygon feature and I can check if it intersects a rectangle:

l = QgsVectorLayer(sys.argv[1], "REGION", "ogr")
r = QgsRectangle(x1,y1,x2,y2)
for feature in l.getFeatures():
    if (feature.geometry().intersects(r)):
        return True


texts: Finding the Point of Intersection Between Two Line Segments and Level of Error

texts:
I am working on calculating the intersection of two lines. Each line segment is represented as two geographic coordinate pairs. Additionally, each line segment is less than  1km in distance.

My question pertains to whether it's feasible to use the standard mathematical calculation for determining the intersection point between two lines (using the formula for the equation of a line while assuming the line segments are on a two-dimensional plane).

Does anyone know what kind error I am dealing with if I were to take that approach?

Obviously, it will also depend on where these line segments are situated as distortion increases as we approach the poles.

I'm open to suggestions for other approaches as well (perhaps converting to a different coordinate system?). I am writing this in Python (not ArcPy).


texts: Conversion between lat/lon and local coordinate system

texts:
I have local coordinate system, this is the red grid on the image. Red grid XY points measured in meters, and I know:

The 0;0 point of the red grid and the lat/lon coordinates of this one point. (the geo system is the green lines on the image measured in lat/lon). This is the one and only point I know values on booth coordinate system (the blue point on the image).
The angle of local coordinate system relative to latitude in degree.

Is it possible the convert values between these two coordinate system with a math formula? If I know the lat/lon values of the point on the green grid, can I get the XY values on the reg grid, and if I know the XY values of any point (and the angle of reg grid in deg) on the red grid in meters, can I got the lat/lon values on the geo coordinate system?

How can I get the transformation formula - if it even possible?

-------------- EDITED ---------------

Thinking how to solve this problem.

Convert lat/lon values of the common point (X Y) and the another (unknown) point from lat/lon to X1 Y1 values. This is the method to get same unit (meters) for different coordinate systems.
Count the difference between X1 - X and Y1 - Y to get vertical and horizontal distances between common point and the original lat/lon values.
Finally I can use formula from this page: https://en.wikipedia.org/wiki/Rotation_matrix but the angle is 90 - (my rotation angle from latitude) in degrees.

Is it right way?

-


texts: QGIS, Multiple lines and wrap on character?

texts:
How to make wrap on characters in Labeling layer settings.
I tried various solutions, but it works.
Does anyone know how to do?


texts: What are AGF Files, how do we put them in PostGIS?

texts:
I came across term which says AGF, after Googling I found its Autodesk Geo Files, But some other say, Atlas Geo Files, Which one of them is right. I also came to know that It is an extension to WKT format supported by PostGIS.

Is there a way where we could make POSTGIS consume the AGF Files.


texts: GRASS v.overlay.or Error when trying to split a polygon with a line

texts:
I am using the Grass plugin in QGIS for Windows (1.8.0).  I am trying to split a polygon with a line feature. I am using the v.overlay.or function within the QGIS Grass plugin.  I have the line layer as input vector map A, and the polygon as input vector B.  When I run the function I get the following error:

 v.overlay ainput=cutlines@CompTest atype=line alayer=1 binput=buildings@CompTest btype=area blayer=1 output=bldg_overlay operator=or

 ￼Operator 'or' is not supported for type line
 Finished with error


texts: ArcMap not allowing attributes of layer to be edited?

texts:
In ArcMap, I am trying to remove some statistical areas manually from a layer of statistical areas. To this end, I am using the Editor Toolbar. But when I try to choose which layer to edit, it won't let me click "ok". What may be causing this? Is there a more simple way to proceed?

The following screenshot is what I see when I try to allow editing of the statistical areas layer. The "ok" button does not turn black (i.e. is not clickable) no matter which layer I choose.


texts: GeoServer WMS Layer not shown in Web page

texts:
I want to show the below link in a web page created asp.net. Please answer me in detail how to do this using OpenLayers latest version 3.7.1.

http://192.168.0.123:8080/geoserver/MilkywayWS/wms?service=WMS&version=1.1.0&request=GetMap&layers=MilkywayWS:ISG2&styles=&bbox=8759192.0,1209982.125,8759198.0,1209993.25&width=414&height=768&srs=EPSG:4326&format=application/openlayers&

I am using the code mentioned below to show the WMS Layer using OpenLayers 3

var layers = [
  new ol.layer.Tile({
    source: new ol.source.MapQuest({
      layer: 'sat' })
    }),
  new ol.layer.Image({
    extent: [8759192.0, 1209982.125, 8759198.0, 1209993.25],
    source: new ol.source.ImageWMS({
      url: '192.168.0.123:8080/geoserver/MilkywayWS/wms';,
      params: {'LAYERS': 'MilkywayWS:ISG2'},
      serverType: 'geoserver' })
    })
  ];
var map = new ol.Map({
  layers: layers,
  target: 'map',
  view: new ol.View({
    center: [8759195.0, 1209988.175],
    zoom: 4
  })
});


texts:
I have the following problem in ArcGIS:
for a species distribution model I need environmental layers and occurrence points.
The problem is my occurrence data points are in ETRS_89_UTM_Zone_32N and some of my layers are in the some, but others are in GCS_Bessel_1841. Obviously its not very easy to transform between these two, I struggled even after using google for some time now.

How do I transform or project everything into one format?

For the occurrence data I need x/y coordinates.


texts: Export only selected and reordered columns to Excel

texts: I am using QGIS 3.12. I want to export my attribute table to an exciting EXCEL table but I want to arrange the columns in a different order and also to select which columns to export and which not. Is it possible?


texts: QGIS Duplicate Labels still appearing despite suggested changes. Do i need to do an expression?

texts:
I have one field called Streets. In Streets there are road beds with suffix as Pkwy, Street, Avenue, Expressway and so on.

I also want to label some of these road beds.

In Layer styling under placement on Repeat i clicked No repeat, and in rendering I clicked Merge connected lines to avoid duplicate labels?

What do i need to do to remove these duplicates and in their place label some other streets from the view below?


texts: Order of WMS layers is reversed in QGIS?

texts:
I have WMS MapService that publicized with ArcGISServer 10.4:

When I connect the WMS in QGIS the order of layers seems to be reverse:

Does anyone has encountered already similar problem?

I selected only the group not separate layers.

The diffrence between the order of the layers I see already when I display REST Url from ArcGIS Server of the same Map Service:

The GetCapabilites from WMS Service shown wrong reverse order of the layers:


texts:
I create various layers programmatically using ArcPy.

I can adjust how the layers are sequenced using the arcpy.mapping.MoveLayer command.

However, I want to have the legend reflect a different order. I have yet to figure out how to do this in Python.

I have tried using legend.autoAdd=True or False. I think I am close in using:

legend=arcpy.mapping.ListLayoutElements(mxd,"LEGEND_ELEMENT","My Legend") [0]

and then legend.items.something.

I have four legend items A, B, C, and D which appear incorrectly as D, C, A, B.

Can anyone provde a link to code examples or better documentation?

This question was originally posed for ArcGIS 10.0 but another question (Moving legend items up and down with arcpy) using ArcGIS 10.2 has been merged into it so answers using any version from and including 10.0 should be acceptable.


texts:
One year ago, we exported data from OpenStreetMap. We added attributes and cut some segments with JOSM.
Now, we want to update those files (.osm) to the current state of OpenStreetMap.
How can we merge our local data and a planetfile for instance ?

We can copy our attributes to the new planet file (using the OSM ID correspondence), but we loose the new segments we created. What can we do merge both ?


texts: From the help document, spline interpolation is supposed to produce a surface runnin through the exact input point. i have about 100 point and run spline interpolation. i also extract the interpolated value from the spline surface at the 100 input points location and compared with the original value. i found that there are some siginificant difference between the interpolated value and original value. Did anybody came across this before? Why the spline surface does not run through the input points. P.S. i used tension spline with 15 as weight and 12 as searching points, cellsize 20m of 3D analyst of ArcGIS10.


texts: How do I reproject one raster to match another's in GRASS

texts:
I have two raster layers, A and B, which have information about the same location but have different projections and datums.

How can I reproject layer A to match the projection and datum of layer B?

Trying to use r.proj yields:

ERROR: Input raster map is outside current region

When I use g.region to set the region to layer B, I get the same error.

When I use g.region to set the region to layer A, then the input and output projections are the same.

Layer A has a projection of:

Projection: Lambert Conformal Conic
N: 5786515.18365295    S: 5349685.18365295   Res:    30
E: -13519426.48582491    W: -13806586.48582491   Res:    30

Layer B has a projection of

Projection: Lambert Conformal Conic
N:    5087310    S:    4800360   Res:    30
E:     606210    W:     440460   Res:    30

The Lambert Conformal Conic is because this is the projection I selected for the Location. Presumably the specific values of the LCC occurred when the raster layers were imported.


texts: Display a WFS service with OpenLayers 3. Define an add custom projection

texts:
I am new in OpenLayers 3. I am working with GeoServer 2.6.2 and I would like to create a web mapping with OpenLayers 3.

I am working with a custom projection from El Salvador (Central America). It is not an EPSG projection. I show you in Proj4 format:

+proj=lcc +lat_1=14.25 +lat_2=13.3166667 +lat_0=13.783333333333333
+lon_0=-89 +x_0=500000.00004 +y_0=295809.17715 +k_0=0.999967040229754
+a=6378206.4 +rf=294.9786982138982 +units=m  +no_defs
+towgs84=0,105.5,197.2,0,0,0,0

OpenLayers library does not recognize this projection. I should define it before calling the WFS Service. I consulted many examples on the web page OpenLayers, but I cannot solve the problem.

Could someone help me?

I tried applying the code @Gabor mentioned, but I can not make it work. This is the JavaScript code that I have created:

proj4.defs('EPSG:1000', "+proj=lcc +lat_1=14.25 +lat_2=13.3166667 +lat_0=13.783333333333333 +lon_0=-89 +x_0=500000.00004 +y_0=295809.17715 +k_0=0.999967040229754 +a=6378206 +rf=294.9786982138982 +units=m  +no_defs +towgs84=0,105.5,197.2,0,0,0,0");

var projSalvador = ol.proj.get('EPSG:1000');
var view = new ol.View({
        projection: "EPSG:3857",
        center: [-9888991.9947, 1533959.4220],
        zoom:12,
        maxZoom:100
});
var vector2Source = new ol.source.ServerVector({
  format: new ol.format.GeoJSON(),
  loader: function(extent, resolution, projection) {
    var url = 'http://localhost/geoserver/workspace/wfs?service=WFS&' +
            'version=1.1.0&request=GetFeature&typename=workspace:Polygon&' +
            'outputFormat=application/json&format_options=callback:loadFeatures' +
            '&srsname=EPSG:6823&bbox=' + extent.join(',') + ',EPSG:3857';
    $.ajax({
      url: url,
      dataType: 'jsonp'
    });
  },
  strategy: ol.loadingstrategy.createTile(new ol.tilegrid.XYZ({
    maxZoom: 19
  })),
  projection: projSalvador,
});

var loadFeatures = function(response) {
  vector2Source.addFeatures(vector2Source.readFeatures(response));
};

var vector2 = new ol.layer.Vector({
  source: vector2Source,
  style: new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'rgba(0, 0, 255, 1.0)',
      width: 2
    })
  })
});

Browser gives me next problem with this code:

"Uncaught SyntaxError: Unexpected token :"  (SyntaxError in url wfs code). OpenLayers can not give me any feature, like I can see in json code [{"type":"FeatureCollection","totalFeatures":0,"features":[{"type":"Feature","id":"Parcela.9415","geometry"....].

I can not find and fix the problem.

Can anyone help me?


texts:
I am using QGIS 3.12.3 on Windows.
I am loading a tif file (CRS EPSG:3395 - WGS 84 / World Mercator - Projected, 1 Band with a palette with 256 colors), and I am trying to just save it "as it is" to another tif file (that is not the final goal obviously).
When I do that, just selecting a file name and leaving all other options unchanged, I end up with a layer with 1 band and a palette of just 2 colors black and white. I can copy/paste the style in the newly created layer to get my colors back, but I am stuck on how to actually save it into a tif file.
I have tried many things, I know I can save it as an image but I don't want that, I want to keep it as a one band/palette GeoTIFF, using the same CRS etc.
The final goal is just to crop the tif, keeping the same palette/resolution/projection/etc.


texts:
I use the plugin regularly and have already used every function. But I have often noticed that the output layers are not always intuitively understandable.
I googled on the topic and found no definition of the result layers. Incidentally, there is no information about the output in the manual.

For example, I am interested in the exact definition of accuracyoutputlayer.tif. the values and their range are not clear to me.  Or after a normal classification, some layer for each class are also created. The meaning of the values is not clear to me too.

Do you have any information?


texts: Importing raster with R from File Geodatabase?

texts:
How can I import a raster file from a gdb into R?
Name of raster file: Cont_maturite
in gdb : CAROPER_04151.gdb
In wd: c:\test\caroper


texts: gdalwarp crop_to_cutline generates error for some states

texts:
I am trying to crop states of a country by using gdalwarp . For some of the states I get this particular error, while for other states it runs fine. This is my second day at GIS/gdalwarp, so please pardon my naiveness. As of now, I do not know anything about shapefiles or geojson (other than how to view them)

Here is the sequence of commands:

#Reprojecting the world population obtained from NASA:
gdalwarp -t_srs EPSG:4326 shapefiles/World/gpw2000.tif inter/repropop2000.tif

#Reprojecting the shapefile containing states of a country (India)
ogr2ogr -t_srs EPSG:4326 India/reproj.shp  India/India_state.shp

#Creating geojson for a particular state
ogr2ogr -f geojson India/state1671.geojson India/reproj.shp -sql "SELECT * FROM reproj WHERE state_id = '1671'"

#Cutting the map
gdalwarp -crop_to_cutline -overwrite -cutline India/state1671.geojson repropop2000.tif lastcut.tif
Warning 1: Too few points in geometry component at or near point 82.551917767195008 19.199543505051185
ERROR 1: Cutline polygon is invalid.

I am using SRS EPSG = 4326. I believe tifs should be fine as the last command works for some other states.

Let me know if I need to attach any of the files.


texts: How to split attribute field of a view in ArcGIS Enterprise Portal

texts:
I am working with a view in ArcGIS Enterprise Portal.  This view contains combined information about the same asset, for example, pipelines.  Some of the information has been combined into a single attribute field, separated by a pipe [ | ].  An example would be pipe diameter: "48 | 54" (for a reducer).

My question is how to split this attribute field by the pipe, for searching and filtering purposes.  I am unable to perform a field calculate on a view. It doesn't necessarily need to be split into separate fields, however if that is the only way, that would be helpful to know.


texts:
I would like to translate from TIFF to ASCII,

import os
os.system("gdal_translate -of AAIGrid -a_srs EPSG:25830 " + 'my.tif' + " " +  'my.asc')

But when I do that the ASCII file is moved:
ASCII result

Driver: AAIGrid/Arc/Info ASCII Grid
Files: my.asc
       my.asc.aux.xml
       my.prj
Size is 76, 70
Coordinate System is:
PROJCS["ETRS89_UTM_zone_30N",
    GEOGCS["GCS_ETRS_1989",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (493800.000000000000000,4181600.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  493800.000, 4181600.000) (  3d 4'13.47"W, 37d46'54.22"N)
Lower Left  (  493800.000, 4174600.000) (  3d 4'13.26"W, 37d43' 7.08"N)
Upper Right (  501400.000, 4181600.000) (  2d59' 2.76"W, 37d46'54.29"N)
Lower Right (  501400.000, 4174600.000) (  2d59' 2.81"W, 37d43' 7.16"N)
Center      (  497600.000, 4178100.000) (  3d 1'38.08"W, 37d45' 0.72"N)
Band 1 Block=76x1 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

Original tif file:

Driver: GTiff/GeoTIFF
Files: my.tif
Size is 76, 70
Coordinate System is:
PROJCS["ETRS89 / UTM zone 30N",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","25830"]]
Origin = (493800.000000000000000,4181600.000000000000000)
Pixel Size = (100.000000000000000,100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  493800.000, 4181600.000) (  3d 4'13.47"W, 37d46'54.22"N)
Lower Left  (  493800.000, 4188600.000) (  3d 4'13.69"W, 37d50'41.35"N)
Upper Right (  501400.000, 4181600.000) (  2d59' 2.76"W, 37d46'54.29"N)
Lower Right (  501400.000, 4188600.000) (  2d59' 2.72"W, 37d50'41.42"N)
Center      (  497600.000, 4185100.000) (  3d 1'38.16"W, 37d48'47.85"N)
Band 1 Block=76x26 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

In ASCII, the file pixel size 'y' is -100 when should be 100. How can I fix this?


texts: Use HTML button tag to select different feature layers in ArcGIS JavaScript API 4.0

texts:
I am trying to build a row of buttons on the top of a page which changes the feature layer being viewed. I am using sample code from here: https://blogs.esri.com/esri/arcgis/2016/01/25/3d-visualization-working-with-objects-paths-and-extrusion/

Here is the HTML:

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Number of Charters by State</title>

  <link rel="stylesheet" type="text/css" href="Map.css">
  <link rel="stylesheet" href="https://js.arcgis.com/4.0beta3/esri/css/main.css">

      <button type="button">About Charters</button>
      <button type="button">Charters by State</button>
      <button type="button">Charters by County</button>
      <button type="button">Charters by Urbanized Areas</button>
      <button type="button">Individual Schools</button>
      <button type="button">Data Access</button>

Line 51 contains the feature layer which I am currently viewing.  What would be the best process to use so that each time a different button is clicked a different feature layer is displayed?

I am not a novice to HTML but I don't practice my coding skills regularly either.


texts:
I'm using ModelBuilder in ArcGIS Desktop 10.6.
I have a model I've created that takes an input .CSV file with coordinates/projection that re-projects the input coordinates into several other coordinate systems.
My issue is the model will only run up to the first "project" process and every process after that becomes invalidated (becomes white and won't run). Once I reconnect the feature layer to the first "project" process it works and runs the rest of the model. Once I save and close the model or validate the model I'm back to step 1 and the model only works to the first "project" process.
Here's the model after I've got it all working:

This is the model after I validate it or reopen it:

After I reconnect the variable below it works again:

I don't want to have to do this every time I run this model.
How do I get the "Project" tool to not break, and stay validated?
I even made variables for input coordinate system, output coordinate system, and geographic transformation, but it still breaks.


texts: How to get the Parameters for a Lambert conformal conic projection CRS in QGIS 2.6

texts:
I am currently having some trouble adding the CRS to a vector file that I want to re project to WGS84/UTM. The shapefile is using Lambert Conformal Conic Projection; ITRF92 horizontal datum and GRS 1980 ellipsoid according to the metadata.

I will like to know how to define the original CRS since in QGIS Lambert CCP/ITRF92 does not appear as part of the available reference systems. The full definition is:

GEOCCS["ITRF92",DATUM["International_Terrestrial_Reference_Frame_1992",SPHEROID["GRS
1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6651"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Geocentric
X",OTHER],AXIS["Geocentric Y",OTHER],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4914"]]

Projection: Lambert conformal conic
Standard parallel: 17.5
standard parallel: 19.5
longitude of the central meridian: -102.000000
latitude of the origin of the projection: +12.000000
false easting: 2 500 000
false northing: 0

The data's location is in Mexico.


texts: Geoserver api: how to determine a layers respective workspace?

texts:
I am using the geoserver api to find information about layers that are available. I am using:

http://localhost/geoserver/rest/layers.json

to get a dynamic list of all layers.
If I then select a layer with:

http://localhost/geoserver/rest/layers/<myLayer>.json

I get information about the layer.

How can I determine the layers respective workspace?


texts: Fill tool creating large output files from sections of a Raster

texts:
(This raster information for one of the clipped rasters)

I have a raster image representing a landscape that is about 1.3Gb. I clipped the raster into 31 sections which vary in size from 10Mb to 40Mb. Whenever I use the Fill tool on a section, the output raster is 1.3Gb (like the whole original raster). This is too much memory to process. Why is the Fill tool creating such a large output file for just a section of the original raster? Is there a way to reduce this size?

Edit: Additional Information

The Rasters were clipping using the methods in this previous question

Clip raster using shapefile into multiple rasters using ArcGIS Desktop?

Edit 2: Python Script changing the mask environment before using Fill (It works, but still provides the 1.3Gb "Fill" Raster)

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace = (r"Workspace\Path")

    for i in range (31):

        arcpy.MakeFeatureLayer_management("shapefile","featurelayer" + str(i),'"FID"='+str(i))

        arcpy.env.mask = ("featurelayer"+str(i))

        arcpy.ClipManagement(“fullextentraster”, “#”, raster +“clip”+str(i)+ “.tif”, “file_name” + str(i), “0”, “ClippingGeometry”)

        outFill = Fill(“rasterclip”+str(i)+“.tif”)
        outFill.save (r“Save\Location\Folder\rasterfill”+str(i)+ “.tif”)

Edit 3: Python Script changing the processing extent environment before using Fill

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace = (r"Workspace\Path")

    for i in range (31):

        arcpy.MakeFeatureLayer_management("shapefile","featurelayer" + str(i),'"FID"='+str(i))

        arcpy.env.extent= ("featurelayer"+str(i))

        arcpy.ClipManagement(“fullextentraster”, “#”, raster +“clip”+str(i)+ “.tif”, “file_name” + str(i), “0”, “ClippingGeometry”)

        outFill = Fill(“rasterclip”+str(i)+“.tif”)
        outFill.save (r“Save\Location\Folder\rasterfill”+str(i)+ “.tif”)


texts:
I downloaded GEBCO bathymetry data and trimmed all my environmental data with ArcGIS according to my GEBCO mask, but Maxent can't run when I insert the GEBCO data. When I delete it and only run it with my environmental layers it seems to work fine. I have already checked that GEBCO data is in .asc form using ArcGIS and this seems to work, but Maxent does not run.

This is the error which appears:

  For input string No data Value