texts: Custom Projection in Proj4 has errors, what can I do to implement a 4 parameter transform?

texts:
I'm working on a project where I'd like to define a projection in QGIS, using a Proj4 string, so that a user can import and export site coordinates from the local UTM or geographic coordinates.

I've had a good search around the place, and the best resource I found is this gis.stackexchange post, which gets me close, but I have errors still, which you can see in the images below. They are not systematic in a way that I can see, but might be related to scale.

This is my Proj4 string:

+proj=omerc +lat_0=-42.83254553 +lonc=147.31794080 +alpha=43.297962 +k_0=0.999592 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

And I have used a bit of software to generate transformation parameters from known coordinate relationships, which results in these:

X shift: 525985.101348
Y shift: 5257731.440819
Z rotation: 0.755692
Scale: 0.999592

These were what I used to get the parameters in the Proj4 string, including coming up with the origin Lon and Lat (via transforming the 0,0 point back to MGA and thence forward to WGS84).

If I use the above parameters to do a manual transformation, the coordinates are transformed perfectly. If I use the above projection in QGIS I'm still getting errors of up to 2 m.

I wonder if anyone knows of a better projection string to use in Proj4, or if there is another way that can handle a site datum or custom coordinate system with a transform into a known system.

EDIT: more info

Here's some coordinate pairs. Note that the lat/lon are derived from the MGA coordinates. E and N mga is GDA94, MGA Zone 55. The Zs are irrelevant it's a simple vertical shift.

Note that I worked out the Lat/Lon shift for the Proj4 string by converting a 0,0 site coordinate to MGA, and used QGIS to convert this to WGS Lat/Lon. This give the origin Lat, Lon as per the other stackexchange post.

Emga,Nmga,Eplant,Nplant,Lon,Lat,Zmga,Zplant
525631.34,5258107.55,-515.61,31.13,147.3135953,-42.82917059,10.7,10.02
525634.87,5258111.13,-515.49,36.17,147.3136383,-42.82913824,10.33,9.65
525633.12,5258109.03,-515.33,33.44,147.313617,-42.82915721,12.43,11.74
525633.98,5258109.93,-515.32,34.68,147.3136275,-42.82914907,9.76,9.08
525692.84,5258084.4,-454.95,56.47,147.3143487,-42.829377,3.09,2.4
525692.41,5258083.45,-454.61,55.49,147.3143435,-42.82938557,3.14,2.46
525856.42,5258016.48,-289.24,119.25,147.3163532,-42.82998312,3.2,2.52
525876.22,5257998.88,-262.762,120.021,147.3165963,-42.83014095,9.59,8.91
525907.25,5258001.29,-241.82,143.07,147.3169758,-42.83011819,3.99,3.3
525906.52,5257997.55,-239.79,139.84,147.3169671,-42.8301519,3.23,2.54
525912.12,5257997.14,-235.43,143.38,147.3170356,-42.8301554,3.17,2.49
525925.34,5257969.74,-207,132.5,147.3171986,-42.8304017,3.41,2.73
525939,5257962.35,-191.98,136.49,147.3173661,-42.83046778,2.89,2.21
525934.89,5257911.72,-160.24,96.81,147.3173181,-42.83092385,2.91,2.23
525933.72,5257861.72,-126.8,59.6,147.3173061,-42.83137416,2.91,2.23
525981.79,5257875.04,-100.93,102.28,147.3178936,-42.83125257,41.97,41.29
526006.26,5257860.79,-73.34,108.7,147.3181937,-42.83138007,8.59,7.9
526086.76,5257842.88,-2.44,150.89,147.3191794,-42.83153861,8.79,8.1
525788.67,5257525.66,-1.84,-284.59,147.3155468,-42.83440535,10.52,9.83
525788.87,5257456.4,45.82,-334.89,147.3155524,-42.83502905,10.47,9.79
526149.54,5257795.97,75.46,159.8,147.3199497,-42.8319589,10.48,9.8
525789.24,5257402.47,83.09,-373.9,147.3155594,-42.83551469,42.49,41.81
526184.03,5257761.77,124.03,158.57,147.3203733,-42.8322657,37.11,36.43


texts: Creating label coordinates in a virtual columns in QGIS 2.18

texts:
I'm working from a database (PostgreSQL/ PostGIS).
I created a view in QGIS of joined polylines from one table with the ST_union function.
I want to label this polylines in QGIS. So, I created 2 virtual columns with $x and $y expressions for the labels coordinates and configured the properties.
But, the virtual columns are not editable and I can't move the labels.

Do you have an idea ?


texts:
I have started having a problem with georeferencing.

When I try to add control points, the cursor just starts "thinking" forever, and it's nigh unusable.  I've read forum posts that say this is an issue with GDI object leaks, and the advice is to reduce the number of rasters in the TOC.  I only have one and it's a real small .jpg, so I just don't understand what's going wrong.

What could have changed on my system to cause such a drastic performance slump?

I am using 10.1 right now, though it was rather recently that I upgraded from 10.0.  I'm fairly certain that I've georeferenced images since the upgrade, but maybe I haven't had to yet.

Is there a significant difference in how georeferencing works in 10.0 compared to 10.1?

The .jpg is stored natively (on the C drive), and actually, I think I did only try one raster.  It was just such a normal little jpeg I didn't think anything of it.  I should try another. I think Snapping may be the answer, actually.  I was trying to snap to a huge point layer, and even though it didn't take long to draw (a misguided proxy for overall performance I realize), I think you're probably right about the issue.

I'll be able to fix this up tonight and give an update.

Yes, it was snapping, the point layer I was hoping to snap to was huge.  Easy fix.


texts: Cache tiles for offline use of Google/Bing aerial maps

texts: I am currently working on a project where QGIS is used for creating and editing vector data and requires a Google or Bing aerial photo layer as a basemap.  I can add Google satellite/Bing aerial photos as a base map layer using the OpenLayers plugin but I am looking for any solution that I can use to do this offline using local cached tiles.


texts:
I'm trying to calculate the logarithm of each cell in a raster. Is there a way to do this with F/OSS GIS? gdal_calc.py only has basic operators and the QGIS raster calculator feature doesn't seem to have a log function either.

I suppose I could still look into PostGIS raster functions but I'd prefer something that I could easily visualize in QGIS.


texts: 'Unknown exception' when editing vector layers in QGIS

texts: I keep coming across the unknown exception error message when trying to edit or add to a vector layer. Anyone else come across this problem? Any fixes? I'm currently running QGIS 3.4


texts: Get LatLng coordinates from SIRGAS2000 UTM 21s coords

texts:
I have a bunch of locations (36k) and they are represented using SIRGAS2000 (UTM 21s) coordinates. I need some way of converting them to traditional LatLng format so I can upload them to my database and use those coords with Google Maps API.

I could use either a formula to translate them to LatLng before uploading, or some kind of API to convert them before using them in my app (after uploading them as they are to my database).

Any clues?


texts: Can I add KML to a georeferenced PDF and retain georeferencing?

texts:
I would like to be able to take a KML file and combine it with a georeferenced PDF so that the KML data shows on top of the georeferenced PDF map within the PDF file.  I have read that i can add .shp files to a PDF but have not managed to accomplish that, but is there a way to not have to convert the KML to shapefile before importing into the PDF?
I need to be able to do this without ArcMap or QGIS.
Thanks.


texts: Reverse geocoding TIGER road data to include US place names

texts: I downloaded primary, secondary, and local road data in Montana via USDA; however, place names (city/town etc) are not included in the data sets. I'm just wondering if there is a method to reverse geocode placenames or if there is road data by state that includes placenames? I looked around for place boundary data but was unable to find anything (I'd spatially join the place name to all roads within the place boundary).


texts: Shapefiles not showing in OpenLayers using QGIS export to OpenLayers 3 plugin?

texts:
I am new to GISs in general and I am working with QGIS to generate layers to be shown on OSM using OpenLayers.

My problem is the following:

I have multiple layers in a QGIS project and I want to export all of them to be displayed in OpenLayers.

To do that I used two methods:

1) Export to OpenLayers 3 plugin of QGIS

The exporting process is done without errors but when I try to visualize some layers they don't appear even when I hide the rest of the layers including the base OSM layer.

2) js-shapefile-to-geojson javascript API

This API converts shape files to geojson format but projecting them on an OSM map doesn't place them correctly.

They appear like a dot in a placement too far from the desired location.
However, when I tried to visualize the same layers that does not appear using the test page provided by the API I was able to visualize them.
So here the problem is to integrate the layers with OSM.

Are there any suggestions on how to visualize these layers using QGIS plugin ?


texts: Reproject MODIS level 3 chlorophyll data using gdal

texts:
I am trying to reproject MODIS Level 3 chlorophyll data (for instance, this file) using gdal in R.

I used gdal_translate and gdalwarp functions. There was no crs information associated with A20121832012213.L3m_MO_CHL_chlor_a_4km.hdf file.

I used:

s_srs= "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +units=m +no_defs +R=6371007.181"

gdalwarp(source,dest,
s_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +units=m +no_defs +R=6371007.181",
t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",srcnodata=-32767,dstnodata=NA)

Although gdal_translate can translate hdf file to other formats, like img, this file cannot be reprojected by gdalwarp because there is an error:

  GeoTransform values not available which is returned by function
  GetGeoTransform

PS:
Description of Level 3 data says:

  This image is typically a scaled value (16bit), two-dimensional array of an Equidistant Cylindrical projection of the globe.

However MODIS Level 3 data is different from SeaWifs Level 3 data. There is no problem with reprojecting SeaWifs data, but not MODIS.


texts: Calculating direction of polygon using other features? (QGIS)

texts:
In QGIS 2 or 3 I have some DEM rasters, on top of that vector layers with plots and buildings + pin (blue one) with coordinates of building.

How to calculate position of garden? Is it n, ne, e, se, s, sw, w, nw? In this case it should return of course: East.

Suppose I can calculate distance from blue pin to all corners of shape (which might be more complicated), choose 2 the longest distances and then somehow calculate direction?

How to do that feature-by-feature? In this case feature means pink plot.


texts: Breaking Mosaic into Large Tiles using GDAL?

texts:
I have 607 aerial images of Massachusetts in .jp2 format.  They are in a mosaic however I need to break down the 607 images into groupings of 16(4x4) spatially so that I can more easily process Massachusetts as a whole.

I've been told that using gdal could help but I am not sure.

Has anyone succeeded in a similar task?


texts:
I have a table of roads and a table of points. I'd like for each point to have the ID of the closest road, as well as the distance to that road.

The line data has both the path type and the geometry type for its coordinates. The points have XY coordinates. What would be the best way to go about finding the closest road using PostGIS?


texts: How to add a few points to an external esri map using iframe and javascript

texts:
There is an external map on : https://gis.apcrda.org/lps/index.html

i have iframed it using this code in my html page :-

      <iframe
        id="crdamap"
        width="500"
        height="400"
        frameborder="0"
        scrolling="no"
        marginheight="0"
        marginwidth="0"
        onLoad="load_plots_for_sale();"
        title="Avenue Plots for Sale" src="https://gis.apcrda.org/lps/index.html">
      </iframe>

the function to be executed on load of iframe should go like this:

function load_plots_for_sale(){
    $(map).on("load", function(){
       console.log('xxxxxxxxxxxxxxxxxxxzzzzzzzzzzz');
    });
}

but i am not able to pick the map, is there any other event which i should listen to? like:

layer.loadend

if so how can i load another layer to the map after all original layers have finished loading.


texts: How to do data augmentation in ArcGIS API

texts:
I got training data for my deep learning model by using export training data for deep learning tool in Arcgis Pro. I would like to do data augmentation on this training data. Could you please suggest the code to do data augmentation? I try to do this by the code below but it has error.
data=arcgis.learn.prepare_data(r'G:\data_training', class_mapping={0: 'tree'}, chip_size=640, val_split_pct=0.1, batch_size=2,transforms=True,resize_to=800)

I had the error:

AssertionError                            Traceback (most recent call last)
In  [61]:
Line 2:     class_mapping={0: 'tree'}, chip_size=640, val_split_pct=0.1, batch_size=2,transforms=True,resize_to=800)

File C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\site-packages\arcgis\learn_data.py, in prepare_data:
Line 1570:  data = (data.transform(transforms, **kwargs_transforms)

File C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\site-packages\fastai\data_block.py, in transform:
Line 504:   assert is_listy(tfms) and len(tfms) == 2, "Please pass a list of two lists of transforms (train and valid)."

AssertionError: Please pass a list of two lists of transforms (train and valid).


texts: What causes large shapefile file size and can you effectively reduce it?

texts:
I have a shapefile that contains all counties of Poland in one shapefile. The shapefile is provided by central cartography agency. The file with counties they provide is around 30MB.

The same shapefile from GADM is ~2MB large and contains the same level of detail (counties) with almost the same attributes. The only information I'm missing from GADM is the actual county identifier attribute that I want to use to match the counties with my data.

I'm running QGIS and I tried to simplify geometries to reduce the size. I've also tried SimpliPy, Generalizer but I can't get more than 0.5% - 1% size reduction. I Googled and searched Stackexchange but didn't find an answer to my question. I also couldn't find help understanding why the size of shapefiles of the same area can vary so much.


texts: Obtaining Connectivity Matrix from ArcGIS Network Analyst?

texts:
Attempting to extract a binary connectivity matrix from network dataset created from river network shapefile.

Any ideas on how to accomplish this?


texts: Make ArcMap label using the first letter from each word in an attribute

texts:
How do I create labels using the first letter of each word from a specific attribute, preferably using Python as the parser?

For example, the attribute name is Common_Name and the values are eastern white pine, balsam fir, plus many others. Cases may be mixed (e.g. some might be Eastern White Pine or Balsam fir), but I'd like the labels to be all capitols, like EWP and BF, respectively.


texts:
I would like to create a dynamic buffer based on contributing drainage area for a stream. I have done some research and I think the best way would be to use the flow accumulation tool but I'm not sure how to isolate just the stream bed pixels of the flow accumulation raster.

I eventually want to create a buffer based on the acreage of the contributing drainage area for the stream so the buffer would get bigger by pre-defined intervals as the drainage area reached certain values.


texts: Getting Exact Latitude and Longitude GMaps.Net

texts:
I am using WinForms alongside with GMaps.Net.

I wanted to get the exact latitude and longitude of the marker so I used the OnClickMarker handler.

I tried using myVar.Location where it only returns the X, Y values according to the scale of the map and not the latitude and longitude of the actual location.


texts: How do I from convert from street level address to block/block group/census track?

texts:
I just downloaded block-level census data, using this answer Where to get 2010 Census Block data?

It output a spreadsheet with a bunch of rows of the format

GEO.id                   | GEO.id2    | Geo.display-label
1000000US22071000100100  | 2.2071E+14 | Block 1002, Block Group 1, Census Tract 1, Orleans Parish, Louisiana

I want to convert "Block 1002, Block Group 1, Census Tract 1" to something like 700 Block of Jefferson St. Etc.

How do I convert these "addresses" into street addresses or longitudes/latitudes?

Edit

I can go individually from street address to census block using the advanced search on American Fact Finder, but I am looking for a way to do a batch conversion.

The FCC has an API that will do block conversions from lat/longitude to blocks: http://www.fcc.gov/developers/census-block-conversions-api


texts: Error while trying to sync - ArcGIS Collector App

texts: One of my staff is trying to sync edits from the ArcGIS Collector App and keeps getting the attached error message. I, and many other staff, are able to sync edits without issues. She has tried several times on WiFi and on LTE (with sync enabled with data) without success. Anyone have any ideas about this error message?


texts: Clipping Census tract shapefiles to boundaries of City Council shapefiles using ArcGIS Desktop?

texts:
I am currently working on a project where I am using Census tract data to figure out the population characteristics of the city council districts of several major US cities. For the cities that I have completed so far, I have been simply uploading the KML files for both the city council districts and the corresponding Census tracts (clipped to the city boundaries when I can find them) to Google Maps. I use Google Earth when the files exceed the 5MB limit on Google Maps. Then, I go through the tracts one by one and in a spreadsheet I record which city council district(s) they are in.

As you can guess, it takes an excruciatingly long time to do this manually.

I was wondering if there was a way to run some sort of script in ArcGIS Desktop that could do this automatically, recording the percent of the area of each Census tract contained within each district?


texts: How to create a polygon with specified dimensions?

texts:
I'm editing a shapefile in QGIS and want to draw in a rectangle with specified dimensions(Say 101ft by 305ft for example).

Is there a way in QGIS to make draw a rectangle with out actually have to use the mouse to manually place each of the points on the map?

What i'm really looking for is ArcMap type functionality where the users places the first point on the map with the mouse and then the user specifies the direction and length of each point relative to the first on.

Thanks!


texts:
I use to transform data points from UTM (WGS84) to WGS84 EPSG:4326. Basically from UTM to Geographic coordinate system both in (WGS84). I tried to do it again but an error message appears.

Translation:

  Failed to export vector file. Error: Error writing object: Error
  creating space object ( OGR error : Failed to create ./A.csv :
  Permission denied )

Original:

  Falló la exportación a archivo vectorial. Error: Errores de escritura
  de objeto espacial: Error de creación de objeto espacial (error de
  OGR: Failed to create ./A.csv: Permission denied)


texts: How can I utilize NumPy arrays to optimize big data geoprocessing?

texts:
I'm interested in learning how to utilize NumPy arrays to optimize geoprocessing.  Much of my work involves "big data", where geoprocessing often takes days to accomplish certain tasks.  Needless to say, I am very interested in optimizing these routines.  ArcGIS 10.1 has a number of NumPy functions that can be accessed via arcpy, including:

NumPyArrayToFeatureClass (arcpy.da)
RasterToNumPyArray (arcpy)
TableToNumPyArray (arcpy.da)

For example purposes, let's say I want to optimize the following processing intensive workflow utilizing NumPy arrays:

The general idea here is that there are a huge number of vector-based points that move through both vector and raster-based operations resulting in a binary integer raster dataset.

How could I incorporate NumPy arrays to optimize this type of workflow?


texts: What are the types of bus routes data, which are made available from Open Street Map (OSM)?

texts: With regards to bus routes, what are the types of data which are made available for download. Is it possible to export a region's bus routes in open street map (OSM) or GIS data format?


texts:
I have received a shapefile with measurement points to be analyzed, but no matter how I try, I can't find how to convert to lat lon. I suppose it would be good to know the projection but when I look into the file (with R and summary()) it says "Is projected: NA" and "proj4string : [NA]".

One measurement-point I know is close to the border Peru-Chile has x,y = 354521, 7997417.8.

How can I find out how to convert this point (and many others) to lat lon?

I tried to follow the instructions of similar questions here, but to no avail.


texts: How to get around the 1000 ObjectIDs limit on ArcGIS Server

texts: I need use the ArcGIS REST API to query my ArcGIS Server services, the where parameter is like ObjectID in (1, 2, 3,...,2000). When the number of object ids is larger than 1000, the returned result shows "Unable to complete operation".


texts: How is a geographic heat map different from a colored statistical thematic map?

texts:
From Wikipedia,

  A heat map is a graphical representation of data where the individual
  values contained in a matrix are represented as colors.

Apparently, a geographic heat map is only one kind, and this is the only kind relevant to me. Cartography is about geographic maps. In cartography, I often distinguish between topographic and thematic maps.

From Wikipedia,

  A thematic map is a type of map or chart especially designed to show a
  particular theme connected with a specific geographic area.

In thematic cartography, we can at least distinguish between types of thematic maps: statistical, color, choropleth, isopleth, dasymetric, point-symbol, and so on.

How is a geographic heat map any different from a colored statistical thematic map?


texts: interacting with GIS without the use of eyes

texts: I've been accepted to do gis, however, I'm not sure what tools are available to blind persons who are in the field. I would like to get information about various methods or tools that blind persons can use to study/work in the area of GIS. There are many ways for a blind person to read maps by touch, i however would like to interact with a map on the level of creating one digitally on the computer.


texts:
I'm using the following SLD for Rendering a Polygon using GeoServer:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
                       xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
                       xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
                       xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>City</Name>
    <UserStyle>
      <Title>City</Title>

      <FeatureTypeStyle>
        <Rule>

          <PolygonSymbolizer>
            <Stroke>
              <CssParameter name="stroke">#FFFFFF</CssParameter>
              <CssParameter name="stroke-width">3</CssParameter>
            </Stroke>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>

      <FeatureTypeStyle>
        <Rule>
          <Title>City</Title>
          <PolygonSymbolizer>
            <Stroke>
              <CssParameter name="stroke">#000000</CssParameter>
              <CssParameter name="stroke-width">1</CssParameter>
            </Stroke>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>

    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

When I make a GetLegendGraphic request, I get an image like:

As you can see, it is showing both the FeatureStyles, because of which the top part of the Legend Graphic looks empty.

I know that there is a FEATURETYPE Parameter on the GetLegendGraphic request, but I can find an example of what I need to specify here.

How do I use the FEATURETYPE for GetLegendGraphic request?


texts: Cannot print the correct values of array in C - GDAL

texts:
I have the following code, which opens a raster file using C, then gets some attributes, and then declares in line :23 a variable called line, then I read the raster as array using GDALRasterIO.

#include "gdal/gdal.h"
#include "gdal/cpl_conv.h"
int main()
{

    GDALDatasetH raster;
    GDALAllRegister();
    raster = GDALOpen( "/home/roger/Documents/98_GDAL_C++_snippets/data/Isle_wight.tif", GA_ReadOnly);
    if(raster==NULL)
    {
    printf("Invalid Raster Dataset");
    }
    else
    {
    GDALRasterBandH raster_band;
    int XSize, YSize;

    raster_band = GDALGetRasterBand( raster, 1);
    GDALGetBlockSize( raster_band, &XSize, &YSize);
    //int *pafScan;
    int nXSize = GDALGetRasterBandXSize (raster_band);
    int nYSize = GDALGetRasterBandYSize (raster_band);
    int *line = CPLMalloc(sizeof(int) * nXSize * 1);
    //pafScan = (int*) CPLMalloc(sizeof(int)*nXSize);
    GDALRasterIO(raster_band, GF_Read,0, 0, nXSize,1, line, nXSize, 1, GDT_Byte, 0, 0);
    // looping through the pixels of this line;
    for (int k = 0; k < nXSize; k++) {
        // this is where you'd do your processing on the pixel
        printf("%d-", line[k]);
        }

        CPLFree(line);

    }
    return 0;
}

However when iterating and wanting to print the values found in the first row of the array, I get either 0, or very long values (I guess because I am printing trash found in memory at this locations.

According to the documentation of gdal, from the snippet I am trying to use, this line should:

There are a few ways to read raster data, but the most common is via the GDALRasterBand::RasterIO()
method. This method will automatically take care of data type conversion, up/down sampling and windowing. The
following code will read the first scanline of data into a similarly sized buffer, converting it to floating point as part of
the operation

Everytime I compile the small program I get the following:

open_array.c: In function ‘main’:
open_array.c:25:5: warning: ignoring return value of ‘GDALRasterIO’, declared with attribute warn_unused_result [-Wunused-result]
     GDALRasterIO(raster_band, GF_Read,0, 0, nXSize,1, line, nXSize, 1, GDT_Byte, 0, 0);
     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

I think that is because I am not assigning the line 25: GDALRasterIO to any variable (array), and then access the values of this variable by subscripting it in the for loop.

Does anyone know how to solve this problem? Should I declare let's say an array of nXSize and then assign it to the GDALRasterIO?

I have a python implementation that makes what I want to accomplish in C:

import numpy as np
import gdal

path = "data/Isle_wight.tif"

ds = gdal.Open(path)
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
print(myarray.shape[0])
#print(columns)
for j in range(myarray.shape[0]-1):
    for i in range(myarray.shape[1]-1):
        print( myarray[j][i])
```[![enter image description here][1]][1]

  [1]: https://i.stack.imgur.com/jUQwI.png


texts: Imposm full database schema import template

texts:
I noticed that the default imposm template for importing OSM data into postgis lacks some types of features (powerlines for example). I was searching the web to find a complete schema import template (in the form of python classes, like the default template) but couldn't find any, does it exist? It would be much easier to have a complete template, then remove what you dont need.

Thanks


texts: Openlayers 5 zoom control has + but zoom out button has "âˆ’" - any idea why?

texts:
Doesn't seem to matter what browser I'm using.  I have the prebuilt ol.js v5.2.0.  As well as the prebuilt ol.css

Here is the code to create the view, it doesn't get any simpler:

    <link rel="stylesheet" href="css/ol.css" type="text/css">

I inspect the elements of the zoom control in the browser debugger, here they are:

   <button title="Zoom in" class="ol-zoom-in" type="button">+</button>
   <button title="Zoom out" class="ol-zoom-out" type="button">âˆ’</button>

Notice the odd characters where it should just say "-".
A quick google of these characters indicate they may be "Unicode 8722 – the MINUS SIGN"

If so, how do I get my browser to display it properly?  Is this my issue or an Openlayers issue?


texts: Seeking advice on Map Maker Pro versus QGIS?

texts:
I am an ecologist in the process of trying to decide whether to purchase Map Maker Pro or utilize the free open-source QGIS platform. Most of our needs are associated with basic map making.

I am a novice in the GIS world.


texts:
I read a few posts here at GIS.SE that points to useful algorithms in CGAL. I'd like to test algorithms in CGAL such as medial axis transformation/polygon skeleton extraction, but the examples on the CGAL website are using text files or eps for input/output.

And I am just wondering what's the best way to use shapefiles and other common GIS formats in CGAL?

For example, can I convert polgon shapefiles to a common data format that CGAL reads/writes?


texts: Minimum wall-thickness of a non convex polygon with holes

texts:
What is the most efficient way to find the minimal wall-thickness (value and location) of a complex, non convex polygon area including holes? See example of a polygon in blue, with the minimum wall-thickness in red, although in this case the location is ambiguous, if the two adjacent lines are parallel.

So far, we have tried:

Subdividing polygon lines and finding minimum point-point lines inside polygon (brute force, not efficient for complex polygons with >10'000 points)
Delaunay triangulation and finding minimal edges inside polygon. Not precise enough, only feasible if combined with subdivision of polygon lines first. Here's example (Nr 3), where the Delaunay triangulation would not find simplex edges in red but would miss out on the minimal wall-thickness in the green box:
Iteratively increasing erosion buffer to find the minimal inset, where the erosion polygon breaks into multiple parts = half of minimal wall-thickness. Problem is to find the location of the minimum wall-thickness with this approach afterwards. Furthermore, erosion doesn't always break into multiple parts and misses "dead ends". Here's an example (Nr 2) that erodes to a line and does give the wrong minimal wall-thickness:
Finding medial axis first, then search for minimum circle on the medial axis that is covering but not overlapping the polygon area. Edit: Problematic are the many "wrong candidates" on the medial axis: Eg. (Nr 1) circle A would be wrong, circle B would denote the correct minimum wall-thickness:


texts:
I have 2 shapefiles in QGIS (one point layer and one Line layer as shown in Figure). I am trying to save these two layers in one GeoJSON. I have seen some of the previous answers such as answer [1] , answer [2] but they are merging the files into one with UNION functions.

My goal is not to merge the features but to save the two layers in one GeoJSON file. For instance, when I import the GeoJSON file to QGIS it should show the point and line layer seperately as show in Figure.
I used Export Features As / GeoJSON but it only save 1 layer.


texts: SQL query for raster point samping within postgis

texts:
I am trying to update my table with the values from a raster stored within the database.

Edit:
I've posted the answer bellow: https://gis.stackexchange.com/a/63655/6711

The query I'm trying to use is Q

update data.points_test_10000
    set sdi_index = (
            select
                st_value(rast,geom)  from
                    rasters.sdi_4326 a,
                    data.points_test_10000 b
                where
                    a.rast && b.geom
                    ) ;

but I'm getting the error:

more than one row returned by a subquery used as an expression

without the set statments the query works as it should


texts: Finding equivalent API-supported data to the 2000 Census

texts:
I'm having some difficulty choosing Census datasets for a project that involve various SES measures at the block-group level. I'll be looking at things like median household income and percentage of people over 25 who have gone to college. These measures are available in the 2000 Census and I can use API to load them directly into R (the programming language I'll be using for the whole project).

What is the equivalent source for these measures that is more recent than 2000 that work with API?

The 2010 Census is formatted differently from the 2000 Census and doesn't seem to offer API functionality for non-demographic data in the SF2. I've considered the ACS 5-year datasets, but I'm worried it won't be as rigorous as the Decennial Census.

Here's the website that has API info for 10-year census.

http://www.census.gov/data/developers/data-sets/decennial-census-data.html

These are the variables for the ACS5. Helpful, but not ideal.
http://api.census.gov/data/2010/acs5/variables.html


texts: ArcMap 10.5 Desktop: Highway shield test symbols not respecting mask layer about it

texts:
How do I stop highway symbol markers from appearing above my mask layer?

The actual highway layer is below the mask layer in the drawing order (so the actual highway line appears below the mask layer, but the text symbol still appears above the mask for some reason).
I tried increasing the layer's label weight and increasing the mask's feature weight but neither resolves the issue.


texts: Reprojecting Tiles Using PixelBender in a Flex Client

texts:
I am pondering the possibility of somehow reprojecting tiles for display in my Flex application using Pixel Bender or something equivalent.  The tile cache in question here is not owned/maintained by me, so I have no control over its output.  I know I could use ArcGISDynamicMapServiceLayer to reproject to the desired coordinate system, but that will be very slow.

Hopefully the following use case explains the need to do so:  I have a Flex app that is using ESRI ArcGIS Online basemaps in the standard Web Mercator projection.  I use a service that returns polygons around and over the polar regions.  Displaying in the Web Mercator projection looks like crap.

So, I guess my questions are:
1.  Does anyone know of any good examples/tutorials...etc showing how this can be done?
2.  Will this overwhelm the processing power of the client browser?


texts: Area outside of bounding box colored in reprojected wms request

texts:
When making a wms request to a geoserver image mosaic layer in which I specify a crs that's different from the layer's native crs (i.e. reprojection occurs) I get a non white background color in the areas outside the bounding box as in the below picture:

The blue color is the same that is set in my sld for value 0. -9999 is the nodata value.

Here is my request. The layer has a native epsg of 4269

https://mygeoserver/geoserver/wms?service=WMS&request=GetMap&bbox=-2300000,-2200000,2700000,850000&srs=EPSG:2163&layers=climate:tmax,climate:states&styles=climate:celsius_web,&format_options=layout:provisional&width=1600&height=800&time=2017-01-14&format=image/png

Here is my style sld:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
  <NamedLayer>
    <Name>celcsius</Name>
    <UserStyle>
      <Name>celsius</Name>
      <Title>Temp C</Title>
      <FeatureTypeStyle>
        <Rule>
          <RasterSymbolizer>
            <Opacity>1.0</Opacity>
            <ColorMap>
              <ColorMapEntry color="#FFFFFF" quantity="-9999" label="" />
              <ColorMapEntry color="#CCFFFF" quantity="-30" label="-30"/>
              <ColorMapEntry color="#66FFFF" quantity="-25" label="-25"/>
              <ColorMapEntry color="#00CCFF" quantity="-20" label="-20"/>
              <ColorMapEntry color="#0099CC" quantity="-15"  label="-15"/>
              <ColorMapEntry color="#006699" quantity="-10"  label="-10"/>
              <ColorMapEntry color="#003366" quantity="-5"  label="-5"/>
              <ColorMapEntry color="#000000" quantity="0"   label="0"/>
              <ColorMapEntry color="#660000" quantity="5"   label="5"/>
              <ColorMapEntry color="#990000" quantity="10"   label="10"/>
              <ColorMapEntry color="#CC0000" quantity="15"   label="15"/>
              <ColorMapEntry color="#FF0000" quantity="20"  label="20"/>
              <ColorMapEntry color="#FF6666" quantity="25"  label="25"/>
              <ColorMapEntry color="#FFCCCC" quantity="30"  label="30"/>
            </ColorMap>
          </RasterSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

How can I remove the blue color (make it transparent or white) outside of the bounding box?

EDIT: This question was closed as a duplicate, but I tried putting &BGCOLOR=0x000000 in my request as mentioned in the linked solution and that only changes the background color inside of the extent from white to black... it doesn't change the light blue color outside of the extent which is what I need. Here's a second image to show this:


texts:
I'd like to show a loading message (Loading...) that shows while the datagrid is being populated.

It would show when the search button is clicked and disappear when the grid results have populated.

I've seen a few possible options including dojox.widget.Standby however, it looks like there is something that is already embedded in the dojox.grid.DataGrid

I can't find any documentation on how to access it or show it in my app. Anybody had any luck inserting a loading message in their app?


texts: Reshape Matrix and perform Zonal Statistics in Matlab

texts:
I have 2 grids: 1) n by m - at 100m 2) n*10 by m *10 - at 10m. These grids cover the same area.

I have to perform zonal statistics at scale of grid no.2. Cell value of grid no.1 is cumulative for the area, so before I'll do my zonal statistics, I need to reshape grid no. 1 in a way that all resulting cells of raster no. 1 would have fractional value (as per above example, say: 1 cell of value 100 would result in 100 cells of value 1).

After resizing my raster with proper values I need to perform zonal statistics on those 2 grids. In GIS enviroment it would not be a problem for me, but I have written in Matlab maybe 10 lines of code, so it may take some time till I get there. So I thought I'd ask if anyone had any experience in such calculation in Matlab.


texts: Python.dll could not be loaded with PyScripter

texts:
I'm trying to configure PyScripter with QGIS. I've installed it using OSGeo4W and I'm following instructions to setup a pyscripter.bat file by copying content of qgis.bat file. I have the following code:

@echo off
call "%~dp0\o4w_env.bat"
call "%OSGEO4W_ROOT%"\apps\grass\grass-6.4.3\etc\env.bat
@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\apps\grass\grass-6.4.3\lib;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
Start "PyScripter" /B "C:\Program Files\PyScripter\PyScripter.exe" --python27 --pythondllpath=C:\OSGeo4W\bin

However when I run this file I get the following error and PyScripter exits after that.

Where am I wrong?


texts: OSGeo4W shell does not accept cs2cs command if typed, only if coppied from a website

texts:
I tried to type in a cs2cs command and OSGeo4W shell did not accept it. When I pasted it from the text on the web it did accept it (see screenshot - command executed is the one pasted and command resulting with an error is the one typed in). suggestions?

]1


texts: POSTGIS get the max length of the polygon and the average width

texts:
I have in postgis a polygon type table.

I need to calculate automatically the maximum length of the polygon:

And the second one is the average width:

I am doubting rigth now if this is possible as although all my polygons will be nearly rectangular in other cases it can be really ambiguous and it is nearly imposible to distinguish between width and height.

Let's see if someone can shed light on this.

Regards,


texts: Sort feature type order in FME KMZ writer

texts:
I have 4 feature types I am writing to a single kmz file using FME.  Currently, these feature types are written as folders in the kmz file and appear in alphabetical order.

Is there any way to sort these folders so a particular folder is always first regardless of the first character in the folder name?  I know there is a way to sort the features within each folder, but I cannot determine how to sort the feature types themselves within the document.

Thanks in advance for your help.


texts:
I was working using drone imagery (regular camera) and there is this alpha band which is preventing me from doing maximum likelihood classification.

I want to remove the alpha band. So far I used raster calculator, separate the image into 3 images for each band (R, G, B) and them using raster merge to join them again, but this is so time-consuming and circuitous.

Is there any simpler way?


texts: Why does ArcGIS Server return a PNG file by default unless I hit the refresh button?

texts:
I'm trying to access a web service that returns a map image in different formats. It first always returns a .png file until I hit the refresh button. It then returns me the map in the correct format.

Do you know what could cause that?

http://appserver/ArcGIS/services/City/MapServer/WMSServer?REQUEST=GetMap&Version=1.1.1&LAYERS=0,2,3,4,5,7,8,9,10,11,12,13&SRS=4326&FORMAT=jpeg&BBOX=43.999861,32.999861,45.000139,34.000139&WIDTH=500&Height=500&Styles=


texts: Add layer info to points, keeps failing in QGIS

texts:
Within QGIS, I am looking to get an extra column in the points file, describing in which layer the point lies.

All points are within 1 layer:
Layer source >SHP file from shapefile.zip.
Validate layer has returned 100%, this layer is used.
Points is generated from a CSV, add geometry to points has been performed and used.

School latitude longitude
1    5.5856567  50.5983856
2    3.7947097  50.497005
3    6.0334594  50.5937665
4    2.8758204  50.4749299
5    2.9487591  50.5603811

Being new to QGIS, I am trying Vector < Data management tools < Join attributes by location. This always returns "error: feature (xyz) has invalid geometry. Please fix or choose ignore". Ignore however leads to an empty dataset.

The layers come from the government, should be fine. The coordinates are generated by Google Maps. Other questions concerning joins address invalid layers and intersections. Layers have been validated here, no doubles exist, all points are within a layer.


texts: ArcObjects - Possible To Add Anything More Than Text To Balloon Callout?

texts:
With code I've scraped together from various posts on this site, I am able to display a simple balloon callout with text in a standard ArcMap session.  Is it possible to add anything more than just text, for example some html or even something like a windows form?

I am aware of HTML popups - but I don't think they will work for my situation.  Without getting to verbose here, the point locations and associated content I want to display don't reside in any feature class.  Instead, I have an event handler wired up in ArcMap listening to an external source that invokes the handler with point location & content.

Any thoughts or suggestions?

IMxDocument mxDoc = (IMxDocument)m_application.Document;
IActiveView activeView = mxDoc.ActiveView as IActiveView;
ITextElement pTextElement;
IElement pElement;
ICallout pCallout;
IFormattedTextSymbol pTextSymbol;
IGraphicsContainer pGraphicsContainer;

IPoint anchorPoint = new Point();
anchorPoint.PutCoords(0, 0);

pTextElement = new TextElementClass();
pElement = pTextElement as IElement;
pTextElement.Text = "HELLO WORLD!!!!!!!!!!!!!!!!";

ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;
ESRI.ArcGIS.Display.IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;

int x, y;
displayTransformation.FromMapPoint(anchorPoint, out x, out y);
IPoint messagePoint = displayTransformation.ToMapPoint(x, y - 100);
pElement.Geometry = messagePoint;

pTextSymbol = new ESRI.ArcGIS.Display.TextSymbol();
pCallout = new BalloonCallout();
pTextSymbol.Background = pCallout as ITextBackground;
pCallout.AnchorPoint = anchorPoint;
pTextElement.Symbol = pTextSymbol;

pGraphicsContainer = activeView as IGraphicsContainer;
pGraphicsContainer.AddElement(pElement, 0);
pElement.Activate(activeView.ScreenDisplay);
activeView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, pElement, null);


texts:
I am looking for a tool like Google Maps Engine Lite or Mapbox where I can draw a simple (buffer) circle around a point with a known diameter.
Are there any quick and easy web mapping tools to draw a single circle and share?
For example, someone on a forums says, "Show me what a 10 mile ban on kite-flying would look like?"


texts: Why did some of my line features in one of the vector layers moved?

texts: I have a map with line features representing roads. Each road has a code which is linked to a CSV file with information on the road. When I opened the file the other day, I realised the CSV files and map doesn't correlate anymore. It seems the codes and information of the line features in QGIS got mixed up and moved around. Is this possible and what can I do to move it back like it was?


texts:
I have geoserver hooked up to an external WFS for around 50 layers.
Using openlayers, if I select them one at a time the all load quite happily.

However I have an option to select all of them at once, on which it locks up geoserver, in that I cannot even access the web-admin (load average on the server is zero though)

Whilst it is probably not entirely practical, I would like to find some way to mitigate this issue. Anyone have any advice on the matter?

(my current train of thought is implement some way to stagger loading all the layers and do them one at a time or in batches, but am unsure if I can query a layer to see if it has loaded)


texts: How to get rectangle with two values in converting polygon to raster?

texts:
I have 112 vector polygon shapefiles. I want to convert every single one to raster layer but instead of getting layer with same shape as vector polygon, I want to get rectangle in which cells that match vector polygon have value 1, and all other cells in rectangle have value 0. And that rectangle for every raster layer should have same size because I want to overlap and sum all 112 raster layers together to get final raster layer in which I will have value gradation among cells between 1 (where is only one of primary raster layers lay) and 112 (where all 112 raster layers overlapping).

Is there a simple way to do that? I don't know how to get rectangle in converting polygon to raster.

I'm using ArcMap 10.0

Here is scheme what I want:

I have many polygons that overlapping but all are inside of boundary of Croatia shape. And I want every polygon have the same extent when converted in raster (red rectangle). At the end I should get final raster layer that look like this schematic:


texts: how to convert grid from Saga gis to have raster in Arcgis?

texts: i am new in Saga Gis, i created grid- LS factor (from Rusle equation for erosis) in Saga gis and now i need to "put" it to Arcgis. How can i do this? I tried a lot of possibilities and none of them worked. Thanks


texts: Where are the attribute summaries in Join Attributes by Location in QGIS 3.2.0?

texts:
I have recently switched from QGIS 2.18.21 to QGIS 3.0, and as of today, to the shiny new QGIS 3.2.0.

I need to do some spatial joining so I was quite surprised to find out that in QGIS3, it seems, it is no longer a possibility to summarise attributes from the joined layer. The take summaries of intersecting features is just gone. As a result I am left with joined layers, where the input layer was a polygon shapefile and join layer point shapefile, but the point attribute data is not transferred to the joined layer, and joined layer attribute table is shown as empty.

Am I missing something here?


texts: Understanding the s2 library (geometry on the sphere, cells and Hilbert curve) for finding points nearby

texts:
I am having trouble understanding how the S2 library works so things in this question may make no sense, but I have a goal I would like to accomplish. Lets say I have many lat/lng point I want to store in a database. The database would be a simple key/value type database (e.g get("MyCoolKey")->returns->"MyCoolValue" ). Later based of a lat/lng point I want to lookup in the database to get all the nearest points (say within 100km).

If I was doing something similar with geohashes I would convert my lat/lng point into a geohash:

48.669, -4.329 -> gbsuv

Then I would store this geohash in the database. I would do the same for many other points. Then to find nearby points I would query for all points in that geohash and all the neighbor:

gbsvh   gbsvj   gbsvn
gbsuu   gbsuv   gbsuy
gbsus   gbsut   gbsuw

I have the following functions I can use to accomplish this goal of mine using s2: https://github.com/mapbox/node-s2/blob/master/API.md .

From what I understand I need to convert the lat/lng points into a cellid? Then I will be able to somehow get nearby points?

Using these available functions how would I then:

Take a lat/lng point and "convert" it
Then based of some lat/lng point find other points in the database


texts:
I am working on a vector map library where you can add plain image and vector layers to a map. The map itself consists of an infinite coordinate system and each layer is a coordinate system of it's own.
However I'm uncertain how layers should be added to the map coordinate system.

The current implementation places the layer's center around given coordinates and scales the layer to a given size:
map.add(new layer(image="url", cx, cy, width, height)) -> Figure1

Another option would be to place the layer's origin (bottom left corner) to the given coordinates:
map.add(new layer(image="url", ox, oy, width, height)) -> Figure2

The third option would be to let the user provide 2 coordinates i.e. bottom left and top right so that width and height can be omitted:
map.add(new layer(image="url", {x, y}, {x, y})) -> Figure3
This Polymaps example uses the third approach

Question:
What would you recommend: Are there any conventions, how layers (image or vector) are usually positioned on a map?


texts: Using OS MasterMap and Transverse Mercator projection in ArcMap?

texts:
I'm trying to use the Transverse Mercator projection in ArcMap, but cannot seem to locate it in the coordinate system list.

I'm very new to GIS and ArcMap.

I'm going to try to use some data from OS Mastermap, which I believe is in the Transverse Mercator projection.

When I look at that projection, it looks very different to the one I get in ArcMap when I select WGS 1984 World Mercator, so I think I'm doing something wrong. I can't see Transverse Mercator in the coordinate system list either.

What I'm hoping for taken from wikipeia, what I have:

(And I think my datum is wrong as is is saying D_WGS_1984 and not OSGB36.)

How do I set up ArcMap to be in the Transverse Mercator projection?


texts: Change point symbol to polygon using Arcgis 10.2.2

texts: I would like to ask how can I convert a point symbol to polygon. I have the capital city (point symbol) and I want to visualize it as a polygon, for better cartographic representation. I'm using ArcGis 10.2.2.


texts: How to position line label right of a line in pyqgis?

texts:
I am using pyqgis to dynamically label layers and show them. However, I could not set the label position to be right of line as shown in the image below.

How could I do it?

Below is the code:

                layer.setCustomProperty("labeling", "pal")
                line_layer.setCustomProperty(
                    "labeling/fieldName", "measurement")

                line_layer.setCustomProperty(
                    "labeling/placement", QgsPalLayerSettings.AboveLine
                )

                line_layer.setCustomProperty("labeling/fontSize", "10")
                line_layer.setCustomProperty("labeling/bufferDraw", True)
                line_layer.setCustomProperty("labeling/enabled", True)
                line_layer.setCustomProperty("labeling/drawLabels", True)


texts: Selecting more than layer with polygon selection or select by location with different CRS in QGIS 2.18

texts: I want to select more than one layer in the same time, select by polygon not rectangle.


texts:
I've just make clean install of Qgis Lisboa under Ubuntu Precise. Latest Qgis version under Windows XP is very unstable (crashes every few minuter while editing shp file).

After new installation there are some problems for me:

I can't see or add independent repositories. So I can't add them under "Repositories" tab. Can I add them anyway? What I need is Table Manager at least. Possibility of install some other plugins could be usefull as well.
I can't find "Join attributes" under fTools menu (it was the same problem after upgrading to Lisboa under Win XP, but deleting fTools directory and reinstalling this plugin solve problem but only under Win XP).

So, I can't work under Windows due to crashes nor under Ubuntu :( Could you help me to make Qgis full functional under Ubuntu?
TIA

Best regards,
Cyfranek


texts: How to create an Atlas based on point attributes?

texts:
Essentially I am wanting to produce an atlas based on a categorical field in a point layer.

i.e. I have a point layer of childcare providers with the categorical field "Provision". I've categorized each feature in this field with "After School Club" , "Breakfast Club" etc, and I now want to produce a set of maps that iterate through each category and show only the points for each. One map of after school clubs, one map of breakfast clubs etc. The extents might be subtly different.

I could do it one by one but it seems like there should be a way to produce an atlas based on the extents of each category? (I do feel i'm missing something obvious :) )

Or alternatively is there a way of automating the creation of a polygon layer and using that as a hidden coverage for the atlas?

EDIT: I've made a little progress with this  - you can use rule based styling to switch on and off features relevant to the current atlas coverage feature.  it actually works fine if all you want to do is show a different set of points. I'm now looking at tying that back to a color scheme and reactive legend.


texts: Is it possible to "re-buffer" river polygon dataset?

texts:
I've river polygon data representing the water quality. This polygon data were derived from polyline dataset according to a buffer distance in the attribute table. Unfortunately I don't have this line dataset. So my "possibly naive" question is - is it possible to "re-buffer" the polygon data set to have polyline data again? Maybe there is a way to draw the centerline of each polygon?
I prefer to work with qgis 2.6 but also have ArcGIS 2.2 available.


texts: Leaflet info box error on same name properties in multiple layers

texts:
Following the 'interactive choropleth map' tutorial, I am able to create my own custom info box in a map with multiple layers.

The properties called are within different layers. However, I have to make sure the properties have different names for each layer. For example using 'Statename' as a property in two layers (called States0 and States1) will result in an error: the layers will not render.

Having a third layer without the property Statename but with another property called Cityname, that is not a problem.

Here is the code used for the infobox

var info = L.control({position:'topleft'});

info.onAdd = function (map) {
    this._div = L.DomUtil.create('div', 'info');
    this.update();
    return this._div;};

info.update = function (pro) {
this._div.innerHTML = '
Information
' +  (pro ?
    'State ' + pro.Statename + '
City ' + pro.Cityname + ''
    : 'hover over map');};

info.addTo(map);

And to use the infobox these three pieces of code for each layer

function infoFeatureStates0(e) {
var layer = e.target;
info.update(layer.feature.properties);}

To empty the infobox

function resetInfoState0(e) {
info.update();}

Eventlistener

function onEachFeatureState0(feature, layer) {
layer.on({
    mouseover: infoFeatureState,
    mouseout: resetInfoState,});}

And similar functions for the other layers: infoFeatureStates1, resetInfoStates1, onEachFeatureStates1 and infoFeatureCities, resetInfoCities, onEachFeatureCity.

Since the info box only shows the property values from the topmost layer visible/accessible (mouseover), I don't see why Leaflet doesn't know which value to return. If States0 is the topmost visible layer, it should return the value for Statename from the feature in Stats0 on mouseover.

It should be easy:
Mouseover to 'select' a feature (which is part of just one layer), then update info box for properties as defined in this._div.innerHTML and return the values if those properties are indeed part of that specific layer.


texts: How to rename points by clockwise direction or change starting point in polygon features?

texts: I have a polygon feature class. I convert it to a point. Name of field (for example !Number!) in a point feature converted from polygon "A2"(The first point should be in the upper left corner (relative to north) of the polygon.) must match the direction of points in a clockwise. I know "AutoIncrement ()" at using Field Calculator, to change the numerical name of points by OBJECTID, but the point ID identifier sometimes do not match with the clockwise direction. Only one method that I know is to copy the contents of a table in Excel and renumber as I need. Then insert into ArcGIS the modified table in Excel and "display XY data"->"Data"->"Export data to point features". But maybe someone knows how to do it easier?


texts: Maybe I have to add new layer with my PNG watermark at XML style, but how? I already tried to change somthing at XML, but don't understand how to add one global layer.


texts:
I am very new to using GIS software.

I have availed CAD drawings of piping network from contractors. I have converted them to .DXF and imported to QGIS. I am using OSM standard map in the background. My CAD layer is not lining up properly with the map layer.

How do I make things work?


texts: Accuracy Assessment in Erdas Imagine 2015

texts:
I'm trying to achieve accuracy assessment in Erdas, but I've faced some difficulties due to my lack of knowledge about it.

I've finished an object based classification in eCognition, but whenever I try to obtain an accuracy assessment, it always give %100 accuracy with TTA mask, so I gave up.

now, I've decided to get that accuracy results from ERDAS because of high reliability to ERDAS among the users.

Workflow:

I've got classified shapefile and the raster file I've used to achieve classification.
in ERDAS, it asks for an .img file which is a thematic raster file it says.
so how can I generate an .img file from my existing .tif and .shp files? or is it even necessary to execute that process?


texts:
I have a feature class with 1 million records. Which is the faster way to copy data to another feature class:
Feature Class to Feature Class or Insert Cursor?


texts: Clipping by geomtery results in Raster extent slightly different to input geometry?

texts:
The raster using arcgis geoprocessing is cliped with a feature, reading it's geometry and as a result, rasters extent slightly differs.
Extent of geom:

str(row.shape.extent)
'570600,537958652 381252,627022596 570698,472917815 381474,211727455 NaN NaN NaN NaN'

Extent of raster (gp and gdal - both shows same):

str(gp.Describe(output).Extent)
'570600,424994 381252,499999738 570699,424994 381474,499999738 NaN NaN NaN NaN'

DS=gdal.Open(output)
DS.GetGeoTransform()
(570600.424994, 1.0, 0.0, 381474.49999973783, 0.0, -1.0)

Having a raster of 1x1 m. resolution that is minor thing...
The problem comes as I produce multiple smaller rasters of the same area to compare and the ones smaller I extend with gdalwarp to the shapefile geometry (the one I clip with). I need rasters of the same extent to be able to compare them with numpy. As not the clipped raster has one column more then the extended one!

As I firstly produce smaller rasters, it is easier to get extent from the row geometry instead of test clipping.

Could enyone explain why clipped raster gets other extent then planned? And how to force gp.clip_management to take the geometry extension?


texts:
I have been asked to create some sort of reporting of OGC services data, and as well some sort of comparison reports to show comparison.  I have not worked much with OGC services, but in review I felt the WMS GetCapabilities would be the easiest to extract.

To start moving on the project, I parsed the XML to a flat file and created a database, easy enough, but as I look I see there must be a better way.  I have several thousand endpoints to review and need to look at the entire request with time.  What tools are available to assist in querying service request, as I learn this product?


texts: Do I need to use ESRI's Create Spatial Type functionality

texts:
Does anyone know if I need to use ESRI's Create Spatial Type functionality when creating a field in a non ESRI / non geodatabase Oracle 11g database?

I have an older customer database in Oracle.  After enabling the datatype ST_GEOMETRY on the database via Oracle's Locator functionality, I added a column to an existing table and populated each row as a point with relevant x,y & ,srid values.

However, ArcGIS v10.2.1 does not even recognize the column as a valid datatype.  It labels it "unknown" in ArcCatalog and returns invalid data type when trying to build a query layer.

Would this indicate I need to use ESRI's Create Spatial Type in order to get ArcGIS to recognize the dataype/column in this non-geodatabase?

If so, can someone explain what it is doing that simply enabling the datatype directly in ORACLE does not?

Any theories as to why ArcGIS will not recognize this ST_GEOMETRY datatype?!


texts:
I want to rename a field of a shapefile in SAGA. I tried:
saga_cmd table_tools 23 -TABLE %dgm_snap_berg% -FIELD ID_1 -NAME Berg_ID

%dgm_snap_berg% - path to the shapefile
ID_1 - field supposed to be renamed
Berg_ID - the new name I want to give
The code runs, but the field is not renamed. What am I doing wrong? This is what SAGA gives out:
Parameters

Table: DGM_Punkte_snapped_Berg
Field: <not set>
Output: <not set>
Name: Berg_ID

Also, I am sure the field ID_1 exists:


texts: Automatically disabling one layer when another layer is turned on

texts: I wanted to find a way to have one layer automatically turn off when another layer was turned on, but was unable to find a way to do this in ArcGIS Online. Is this operation possible in ArcGIS Online?


texts: Installing spatialite-gui using ubuntugis

texts:
I am unable to install spatialite-gui using the ubuntugis unstable ppa. There appears to be a conflict between the GEOs libraries required by QGIS, which I installed first successfully (GEOS 3.5.0 automatically installed as a dependency), and the GEOS libraries required by spatialite-gui. Here is the error I receive when trying to install spatialite-gui:

The following packages have unmet dependencies:
 spatialite-gui : Depends: libgeos-c1 (>= 3.4.2) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

When I try to install libgeos-c1, I see that this install would remove my QGIS install. Is is possible to have both QGIS and spatialite-gui installed from ubuntugis ppa, either w/ GEOS 3.5 or another version?

Edit:
Ubuntu 14.04.4
QGIS 2.16.0


texts: ArcMap adds decimal to string conversion

texts:
ArcMap dropped the leading zero from the census block ID for California. When I added a new text field, GEOID, in the field calculator with the formula,
GEOID = "0" + str(!BLOCK_2020!)

ArcMap created the text field but populated it with one decimal place -- for example, 060014001001001.0.
Any idea what I'm doing wrong?


texts: OpenLayers feature styling change by the baselayer

texts:
How can I have the style of the features on a vector layer depending on the OpenLayers' "base layer"?

For example, some base layers are dark images and I need (for the contrast) the features to be displayed bright in the vector layer upon the base layer. And as soon as base layer changes to a bright one, features should be redrawn dark.

How can I achieve this with OpenLayers code?


texts: Expression-based labels in scripted Atlas fail to update

texts:
I try to follow Programmatically load composer from template and generate atlas using pyQgis to automate atlas print composer using pyqgis :

specific in atlas code say :

    myAtlasMap = myComposition.getComposerMapById(0)
    myAtlas = QgsAtlasComposition(myComposition)
    myAtlas.setCoverageLayer(atlas_desktop)
    myAtlas.setComposerMap(myAtlasMap)
    myAtlas.setFixedScale(True)
    myAtlas.fixedScale()
    myAtlas.setHideCoverage(False)
    myAtlas.setFilterFeatures(True)

this code work for my task and i take new 100 atlas feutures maps but in my template.qpt i have some labels with code like this :

QGIS [%concat( "OBJECTID" )%]

that OBJECTID is the same for all new 100 atlas feutures maps.
Any idea why? If I use manual atlas export, then for any map OBJECTID change.

I have more details the true Does not work anything related to layer atlas in my template.qpt i have some base styles rules or base labels rules like this :

intersects($atlasgeometry, $geometry )

or

$id = $atlasfeatureid

that is the some in 100 maps :

i think so that Programmatically load composer from template and generate atlas using pyQgis just export one by one my features from my atlas layer and the atlas in my template stay False.

any idea ?


texts: Has anyone found any informative documentation about Esri's ArcPy Python module?


texts: I want to know how can I set a "selector layer" (L.Control.Layers) outside the container in Leaflet. Like this screenshot of an OpenLayers application:


texts: Is PostGIS ST_Distance equivalent to QGIS distance to hub?

texts: I'd like to make the equivalent function of QGIS's distance to hub algorithms. Is ST_Distance the right approach to start with? Or is that algorithm transposable into something I can let Postgres PostGIS?


texts: Splitting polyline by points using ArcGIS 10.2 for Desktop?

texts:
I have a layer containing a curved polyline and I would like to split it at multiple point locations, which I have in another layer. In ArcGIS 10.2, I can do this with the tool "Split Line at Point". It is working, but now I have a chaotic attribute table. The order of each feature in the table is not the same as it is geographically.

How can I order the polylines geographically?

The line is curved at many sections, so order only by latitude or longitude is not possible. (I think that this method is applied by ArcGIS, because the order of the points for splitting is the correct.)
I have tried it with searching for the same End_X and Start_X, but it seems that they have not always the same coordinates.


texts:
I want to do tasseled cap transformation on rapid eye data. the main question would be about the parameters which are used. Somebody has a suggestion?

A short update concerning this attempt to calculate Tasseled Cap transformation.
After investigating in this area from my point of view, the result of a TC calculation for Rapid Eye will be not sufficent, because of simple reasons. Calculating the 6 TC image (Brightness, Wetness, Greenness, Fourth, Fifth and Sixth) is in Landsat obviously based on 6 Bands (1,2,3,4,5,7). It is possible to use the thermal band but most publications do not consider it. Rapid Eyes "highest" band in that way is NIR which is not really the decisive factor. To calculate Rapid Eyes TC are definitively needed 2 more bands above the red-edge band of Rapid Eye. It possible but i believe it would not make so much sense.


texts: Clipping Vector by Custom Raster Extent in QGIS?

texts:
I have a raster image (elevation map from DEM) where I have limited the rendered elevation to only 3,500' and higher.  I would like to convert the boundary/extent of this rendered image into a vector file, so that I can upload as a shapefile into other programs.  I have not been able to do so yet after much Googling and experimentation.

Among other things, I have attempted to clip a large vector by the raster (I get an error saying "Could not load source layer for OVERLAY.").  I have attempted to polygonize a two color raster (black and white) and keep only one of those.  I believe they are all on the same CRS.  I can convert it into a vector, but that converts all the contours, which is not my goal.  I only want the boundaries.

Picture below of the raster extent (in white) that I want in vector form.

How do I clip out the black portion?


texts:
I am currently having slight problem on displaying marker in the leaflet. I have used the following codes on displaying marker and map, however the map is displaying flawlessly but the marker is invisible. I have also tried worldCopyJump true but no help.

var map = new L.Map('map', {
    center: [0, 0],
    zoom: 3,
    layers: [new    L.tileLayer('http://localhost:8080/1.0.0/osm/{z}/{x}/{y}.png',{tms: true})]
});

var marker = L.marker([85.29009, 27.69545]).addTo(map);

//following also won't work
//var marker = L.marker([85.29009, 27.69545]);
//map.addLayer(marker);


texts: When generating a line using the Douglas-Peucker Algorithm the resulting object is a line with steep angles, I would like to "smooth" that line and making it look less "squarish" using curves, like bezier curves: there is an algorithm for that ?


texts: Changing labelling of OpenStreetMap Base Map for ArcGIS Online?

texts:
I am creating a Story Map and using Open Street Maps as the base map.
One issue my superiors see is the labeling of the streets. It's way to frequent!

See the attached screenshot:

Is there a way I can change the frequency of the labeling for the same road?


texts:
I am tring to intersect lines with a polygon grid, but i got an error when i run this:

create table "st11_01_vd" as
  SELECT (
    "grid_1km.gid",
    "st11_012017_ln.linetime",
    ST_INTERSECTION("grid_1km.geom", "st11_012017_ln.geom")
  FROM "grid_1km",
    "st11_012017_ln"
  WHERE ST_Intersects("grid_1km.geom", "st11_012017_ln.geom");

Where I am wrong?


texts: How to extract unique values from .pbf file with GDAL

texts:
I would like GDAL to return the unique values in a column. Reproducible example:
wget https://github.com/ropensci/osmextract/raw/master/inst/its-example.osm.pbf

ogrinfo its-example.osm.pbf
ogrinfo its-example.osm.pbf lines > lines.txt
head lines.txt

## --2021-04-11 09:32:38--  https://github.com/ropensci/osmextract/raw/master/inst/its-example.osm.pbf
## Resolving github.com (github.com)... 140.82.121.3
## Connecting to github.com (github.com)|140.82.121.3|:443... connected.
## HTTP request sent, awaiting response... 302 Found
## Location: https://raw.githubusercontent.com/ropensci/osmextract/master/inst/its-example.osm.pbf [following]
## --2021-04-11 09:32:39--  https://raw.githubusercontent.com/ropensci/osmextract/master/inst/its-example.osm.pbf
## Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.109.133, 185.199.108.133, 185.199.111.133, ...
## Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.109.133|:443... connected.
## HTTP request sent, awaiting response... 200 OK
## Length: 40792 (40K) [application/octet-stream]
## Saving to: ‘its-example.osm.pbf.4’
##
##      0K .......... .......... .......... .........            100% 3.04M=0.01s
##
## 2021-04-11 09:32:39 (3.04 MB/s) - ‘its-example.osm.pbf.4’ saved [40792/40792]
##
## INFO: Open of `its-example.osm.pbf'
##       using driver `OSM' successful.
## 1: points (Point)
## 2: lines (Line String)
## 3: multilinestrings (Multi Line String)
## 4: multipolygons (Multi Polygon)
## 5: other_relations (Geometry Collection)
## INFO: Open of `its-example.osm.pbf'
##       using driver `OSM' successful.
##
## Layer name: lines
## Geometry: Line String
## Feature Count: -1
## Extent: (-1.561196, 53.806303) - (-1.549845, 53.809293)
## Layer SRS WKT:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",

ogrinfo \
  -dialect sqlite -sql "SELECT highway, COUNT(*) from lines" \
  its-example.osm.pbf

## Had to open data source read-only.
## INFO: Open of `its-example.osm.pbf'
##       using driver `OSM' successful.
##
## Layer name: SELECT
## Geometry: None
## Feature Count: 1
## Layer SRS WKT:
## (unknown)
## highway: String (0.0)
## COUNT(*): Integer (0.0)
## OGRFeature(SELECT):0
##   highway (String) = footway
##   COUNT(*) (Integer) = 189

ogrinfo \
  -sql "SELECT DISTINCT highway from lines" \
  its-example.osm.pbf

## Had to open data source read-only.
## INFO: Open of `its-example.osm.pbf'
##       using driver `OSM' successful.
##
## Layer name: lines
## Geometry: None
## Feature Count: 13
## Layer SRS WKT:
## (unknown)
## highway: String (0.0)
## OGRFeature(lines):0
##   highway (String) = footway
##
## OGRFeature(lines):1
##   highway (String) = cycleway
##
## OGRFeature(lines):2
##   highway (String) = service
##
## OGRFeature(lines):3
##   highway (String) = residential
##
## OGRFeature(lines):4
##   highway (String) = trunk
##
## OGRFeature(lines):5
##   highway (String) = steps
##
## OGRFeature(lines):6
##   highway (String) = unclassified
##
## OGRFeature(lines):7
##   highway (String) = track
##
## OGRFeature(lines):8
##   highway (String) = pedestrian
##
## OGRFeature(lines):9
##   highway (String) = trunk_link
##
## OGRFeature(lines):10
##   highway (String) = tertiary
##
## OGRFeature(lines):11
##   highway (String) = (null)
##
## OGRFeature(lines):12
##   highway (String) = corridor

# fails
ogrinfo \
  -sql "SELECT DISTINCT other_tags from lines" \
  its-example.osm.pbf > test.txt # works

# fails
ogrinfo \
  -sql "SELECT DISTINCT hstore_get_value(other_tags, 'bicycle') from lines" \
  its-example.osm.pbf

## ERROR 1: Invalid index : -1

What is the error message and how can I stop it?


texts: Base-R legend for quick plot of sf-objects

texts:
I have two sf-objects that are polygons and one that is a layer of points. The attributes (columns) don't really matter for quick visualization purposes. I just want to make a basic legend and give them different filling colors and maybe a custom name. Something  like:

black   -- "houses"
green   -- "inner city"
red     -- "outer city"

My code example:

library(sf)

# make multipoints
mp_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
mp = st_multipoint(mp_matrix)
mp_sfc = st_sfc(mp)

# make outer poly
bbox = st_bbox(mp_sfc)
outer_poly = st_make_grid(bbox, n = 1)

# make inner poly
polygon = rbind(c(2,2.5), c(4,2.5), c(4, 3.5), c(2, 3.5), c(2,2.5))
poly_list = list(polygon)
poly = st_polygon(poly_list)
inner_poly = st_sfc(poly)

# plot all
plot(outer_poly, border = "red")
plot(inner_poly, add = T, border = "green")
plot(mp_sfc, add = T, col = "black")

produces the following plot. And all I want is just a really basic legend to the side with the colors and a self chosen name. I think that'd be really helpful for me for some future quick plots when ggplot, tmap or mapview is just too much. Maybe it's even easier with them and I just dont know.


texts: Export map with legend using OpenLayers 6

texts:
I would like to export a map with a legend in the same PDF using OpenLayers 6.
I searched in the OpenLayers documentation but got a map without a legend.
How can I do this?


texts: In CartoDB, can columns (e.g. ARPU2011) have descriptions (e.g. Average Revenue Per User 2011)

texts:
I have an in-house database with spatial data (Postgres + PostGIS) for which I'm generating an extract (GeoJSON/KML etc) in order to load it into CartoDB.

My column names are technical (e.g. ARPU2011, or even less clear, and I need this to remain that way) but I also have column descriptions available (e.g. "Average Revenue Per User 2011").

Using which export format and what type of formatting can I then later find both the technical name and the description usable for filtering in CartoDB?


texts: Alternatives to bilinear interpolation for raster extraction in R

texts: I have a set of coordinates for which I want to extract temperature values from a raster. I am using the extract function from the raster library in R. Looking at the interpolation methods, there are options for simple and bilinear interpolation. I am not very familiar with raster extraction literature but I was wondering if there are other alternatives in R for other interpolation methods which give a better estimate for values at a point than bilinear interpolation. Also is bilinear the usual method used in practice? Thanks.


texts:
I have to use 64bit python. So, I install Background Processing 64bit and change the interpreter of Spyder into python 64bit in ArcGISx6410.3.

but I can import successfully in ArcMap' python window

my OS is window7 64 bit
and ArcGIS is 10.3.0


texts: Pyqgis action when something is selected

texts:
I want to make a pyqgis module that will perform an action similar to Identify tool.
For example there is a tiny qt widget with a text box that is opened until user decide to close it manually. When user selects feature(s) of some layer, script immediately reads feature's attributes and print it in a text box. As soon as user selects another feature, widget's text will also change according to current selection.
Is there some onSelectionChange() method for such actions?


texts: How to resample the raster cells to make a merge (mosaic)?

texts:
I am creating a DEM with raster files, but I would like to downsample the pixel values (resampling to lower resolution/larger cell sizes) in order to reduce the size of my generated file. My idea is to do this using the least amount of memory possible. My code is working fine for mosaic, without the downsampling. And I found this documentation for Resampling, but I didn't figure out how to insert the code sample in the mosaic code.
import os
import rasterio
from rasterio.merge import merge
from rasterio.plot import show
from rasterio.enums import Resampling

layers = ['A', 'B', 'C', 'D', 'E']
for layer in layers:
  # files to merge into the DEM.
  files = ['example1.zip', 'example2.zip', 'example3.zip', 'example4.zip'] # In my real problem, each layer has different files.
  files_path = '/content/raster'

  src_files_to_mosaic = []
  if not os.path.exists('/content/'+layer+'.tif'):
    for grade in list(set(files)):
      q = os.path.join(files_path, grade) #path of each grid
      raster_path = 'zip://'+q+'!'+grade[:8]+'.tif'# extract zip
      print(raster_path)
      src = rasterio.open(raster_path) # open each raster
      src_files_to_mosaic.append(src) #joins the paths to do the mosaic

    mosaic, out_trans = merge(src_files_to_mosaic)
    out_meta = src.meta.copy() # Copy the metadata
    # Update the metadata
    out_meta.update({"driver": "GTiff",
                    "height": mosaic.shape[1],
                    "width": mosaic.shape[2],
                    "transform": out_trans,
                    "crs": "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
                    }
                    )
    # Plot the result
    show(mosaic, cmap='terrain')
    # Write the mosaic raster to disk
    raster_name = '/content/'+layer+'.tif'
    with rasterio.open(raster_name, "w", **out_meta) as dest:
      dest.write(mosaic)


texts: Finding Perpendicular Distance and Minimum Perpendicular Vector between Point and Line using GeoTools and JTS?

texts:
I have a line formed by 2 lat/lon pairs, and the lat/lon of a point. I would like to find out the perpendicular distance between the line and the point on the Earth surface (can assume Earth as great sphere), and the minimum perpendicular vector (i.e. the projected "cross point" on the line).

I am trying to use Geotools 8.0 and JTS for this. Below captured my testing code:

    //Coordinates in lon, lat
    Coordinate linePt1 = new Coordinate(-5.71472, 50.06639);
    Coordinate linePt2 = new Coordinate(-3.07000, 58.64389);

    //Multiply all longitudes by the cosine of latitude.
    //http://gis.stackexchange.com/a/29713/10772
    linePt1.x = linePt1.x * Math.cos(linePt1.y);
    linePt2.x = linePt2.x * Math.cos(linePt2.y);

    LineString line = createLine(new Coordinate[]{linePt1, linePt2});

    Coordinate pt1 = new Coordinate(-6, 54);
    pt1.x = pt1.x * Math.cos(pt1.y);
    Point point = createPoint(pt1.x, pt1.y);

    double distanceOp = DistanceOp.distance(line, point);
    System.out.println("Distance = " + distanceOp);

    //Find the minimum perpendicular vector using "closestPoints()"
    for (Coordinate c : DistanceOp.closestPoints(line, point)) {
        System.out.println("=== " + c);
        //verify if the point is on the line
        System.out.println(CGAlgorithms.isOnLine(c, new Coordinate[]{linePt1, linePt2}));
    }

the createPoint() and createLine() methods:

public static LineString createLine(Coordinate[] coordinates){
    GeometryFactory factory = new GeometryFactory(new PrecisionModel(
            PrecisionModel.FLOATING), WGS84_SRID);
    LineString line = (LineString) factory.createLineString(coordinates);
    return line;
}

public static Point createPoint(double longitude, double latitude) {
    if (longitude < -180 || longitude > 180) {
        throw new IllegalArgumentException(
                "Longitude should be between -180 and 180");
    }
    if (latitude < -90 || latitude > 90) {
        throw new IllegalArgumentException(
                "latitude should be between -90 and 90");
    }
    GeometryFactory factory = new GeometryFactory(new PrecisionModel(
            PrecisionModel.FLOATING), WGS84_SRID);
    Point point = (Point) factory.createPoint(new Coordinate(longitude,
            latitude));
    return point;
}

However, the result of "isOnLine()" return false. I was wondering if there is anything wrong.

Is there something wrong with my verification method, or actually the way I used to find out the perpendicular distance and the "cross point" on the Earth surface is not correct?


texts: Extracting an encoded linestring from a website

texts:
I need to extract this highway route on this webpage.
https://scenicbyways.info/byway/80452.html
web scraping is not part of my skill repertoire. Sometimes I add .geojson to the end of map link and it returns the geometry and attributes which i then scrape with python. however that is not the case here. I DO NOT want to manually create this.

https://scenicbyways.info/byway/80452.html


texts: In GeoServer 2.16 or 2.15 colours raster with 5 bands shown as black and white raster

texts: I have problem with 5 bands colour raster. When I open the raster in ESRI ArcGIS or QGIS the raster is open with colours, but when I uploaded it to  GeoServer software the WMS shows it only with black and white colours. I guess I have to do some settings, but I don't know what and where.


texts: I have tried to fill a gap in a polygon in QGIS, but the tool that uses QGIS instead of that draws another polygon over the gap that I want to eliminate.


texts: Getting working hyperlinks in the feature attributes of a WMS in Qgis

texts:
I have a Geoserver WMS of point locations with urls in one of its attribute fields. If I do a getFeatureInfo request on this in Qgis, using identify, it just returns a string without hyperlink functionality on the urls. For my application it should be able to link to external websites and open local file folders. Preferably without having to copy and paste.
So far, I tried to convert this string to a html link by wrapping it in html-tags, like:

open

As shown in the image the hyperlink tries to open in another instance of QGIS instead of a browser, which fails.

My Question

Is there any way to disclose urls in the feature attributes of a WMS in Qgis with working hyperlinks?


texts: Good advice regarding Mobile GIS Hardware and software for outdoor tasks?

texts:
Could anyone offer advice regarding outdoor mobile hardware and software for GPS-photos, data capture, and textual input. We don´t want to necessarily write long text, so the ability to write into shapefile format would be optimal. Obviously the price is a big factor and we have a budget of around 2-3,000 Euro/Dollar for everything. The colleague who needs this also has an iphone...so this could also be an option.

I know that this has been covered in some extent in other posts, but I also would like to know if anyone has used pure mobile GIS opensource software on Android or Apple based operating systems. The Ipad also has GPS build in and would also be a great platform for GIS mobile tasks.


texts: How to calculate proportion of a line segment around a point?

texts:
I am looking at certain variables surrounding points that represent Sea otter raft sites. The coast line surrounding the points is broken into unit lines (Essentially a polyline with attributes for various segments). For example, The unit lines have Oil Residency Index (ORI)(from 1 to 5) as well as a physical exposure proxy for each segment (segments range in length).

I want to calculate the proportion of different ORI and Exposure classes around each point (500m buffer). IE- What proportion of ORI class 1,2,3,4,5 are there around each point respectively.

I am using ArcGIS 10.2

Thoughts and ideas?


texts:
I looked at the other similar questions posed and none of them seem to be answered.

Basically I have a NetCDF file with multiple subdatasets and from what I've read, I should be converting it to a TIFF file, which I can then import into PostgreSQL using raster2psql.  I'm trying to use gdal_translate command but I don't understand how to batch translate all of the subdatasets in one go.


texts: Implementing proportional stratified random sample in R

texts:
I've completed a landcover classification on an Landsat Image and would like to start an accuracy assessment.  I'd like a stratified random sample that has a sample number that varies based on each landcover values total image pixels.  In other words a proportional stratified random sample.

My code so far is the following

spdf = sampleStratified(FinalClassImg, 200, sp = T)

This generates 200 points per landcover class but what I really need is a variable number of samples based on the overall landcover pixel totals.

I roughly calculated the percentages in the image and decided to use 5000 total points.  I multiplied the percentages * 5000 and got the following numbers for the nine landcover classes

100, 439, 643, 172, 194, 45, 39, 27, 3341

These numbers represent the number of random samples I'd like selected within each landcover class.

I would like to feed this set of numbers into the above code in place of the 200

I tried this:

spdf = sampleStratified(FinalClassImg, c(100, 439, 643, 172, 194, 45, 39, 27, 3341), sp = T)

but I get a mess of warnings

1: In if (nrow(y) < size) { ... :
  the condition has length > 1 and only the first element will be used
- fewer samples than requested found for strata: 2, 3, 4, 5, 9

and only the first number is used (100 samples per class).

Is there a way to do this or does the sampleStratified only work with a single number?


texts: Print composer changes maps when I zoom out?

texts:
I am making a small multiples visualization of about 32 different maps, all from the same shapefile/source using QGIS 3.0.1. I am using the "Provide feature filter" to select the data that I need to display.

I update each map, one at a time, but when I zoom out using my mouse, it completely erases all the work I've done and changed every single map to my most recent filter selection.

Prior to zoom:

This is what it all changes to as soon as I zoom out:

In some of the maps, as soon as I zoom out, they lose the ability to "lock layers."

I tried "lock layers" as well as "lock layers" and "lock styles for layers." I've tried using the same CRS for every map as well as selecting "draw map canvas items."

I am at a complete loss.


texts: Removing border that appears in map exported with arcpy.mapping.ExportToPDF?

texts:
I am exporting the current mxd map with all the layers that appear but it saves it with a framework around it. How can I get rid of it? It appears as well when I export to other formats.

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
df.scale = 1250000
arcpy.mapping.ExportToPDF(mxd, "C:/Users/DickVater/Documents/ArcGIS/IASTE internship/Final map/year/" + name_save + ".pdf")


texts: Why the shape of MOD02 image that was processed by gdalwarp is different with the image by the MCTK

texts:
I am trying to project MOD021KM from Sinusoidal to WGS84 using gdalwarp with Python. The original image with HDF format was downloaded from NASA EOSDIS Reverb website. I used the code below. Moreover, I applied the MCTK in ENVI to project the image. The shape of two processed image is different.

gdalwarp -of GTIFF -r near -t_srs EPSG:4326 -ot Float32 HDF4_EOS:EOS_SWATH:"MOD021KM.A2017002.2050.005.2017003201059.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB 2017002_F32.tif

0...10...20...30...40...50Warning 1: Metadata exceeding 32000 bytes cannot be written into GeoTIFF. Transferred to PAM instead.

The result is shown in the figure.

Why the shape of two image is different? How do I make the image which is processed by gdalwarp same shape as the image by the MCTK.

UPDATE:
The gdalinfo is following:

The image by gdalwarp:

Driver: GTiff/GeoTIFF
Files: 2017002_F32.tif
   2017002_F32.tif.aux.xml
Size is 4478, 1647
Coordinate System is:
GEOGCS["WGS 84",
  DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
          AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0],
  UNIT["degree",0.0174532925199433],
  AUTHORITY["EPSG","4326"]]
Origin = (-161.678555808543310,78.155197825739222)
Pixel Size = (0.013771295523432,-0.013771295523432)
Corner Coordinates:
Upper Left  (-161.6785558,  78.1551978) (161d40'42.80"W, 78d 9'18.71"N)
Lower Left  (-161.6785558,  55.4738741) (161d40'42.80"W, 55d28'25.95"N)
Upper Right (-100.0106945,  78.1551978) (100d 0'38.50"W, 78d 9'18.71"N)
Lower Right (-100.0106945,  55.4738741) (100d 0'38.50"W, 55d28'25.95"N)
Center      (-130.8446251,  66.8145360) (130d50'40.65"W, 66d48'52.33"N)
Band 1 Block=4478x1 Type=Float32, ColorInterp=Gray
Description = Earth View 250M Aggregated 1km Reflective Solar Bands Scaled
Integers
NoData Value=65535
Unit Type: none

The image by the MCTK:

Driver: GTiff/GeoTIFF
Files: 2017002_mctk.tif
Size is 2973, 2622
Coordinate System is:
GEOGCS["WGS 84",
  DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0],
  UNIT["degree",0.0174532925199433],
  AUTHORITY["EPSG","4326"]]
Origin = (-167.328659057617190,78.277236938476562)
Pixel Size = (0.023494699706713,-0.008966192967982)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
  Upper Left  (-167.3286591,  78.2772369) (167d19'43.17"W, 78d16'38.05"N)
  Lower Left  (-167.3286591,  54.7678790) (167d19'43.17"W, 54d46' 4.36"N)
  Upper Right ( -97.4789168,  78.2772369) ( 97d28'44.10"W, 78d16'38.05"N)
  Lower Right ( -97.4789168,  54.7678790) ( 97d28'44.10"W, 54d46' 4.36"N)
  Center      (-132.4037879,  66.5225580) (132d24'13.64"W, 66d31'21.21"N)
  Band 1 Block=2973x1 Type=Byte, ColorInterp=Gray


texts: Generating TIN from elevation data in Java

texts:
I would like to generate TIN from standard elevation data in Java. Are there any libraries I can use? I could not find any just by googling. In addition library should be suitable for commercial use and has to have no native dependencies.

Other option is me writing a library that does the job. For that, can you suggest a good source that explains the algorithm clearly.


texts: Error: 'Raster map not found' in QGIS 2.18.1 GRASS tools

texts:
I get the same error in multiple QGIS GRASS tools and I don't understand the error. I have a raster surface height map with some no data spaces I want to fill.

The error occurs when using r.fillnulls, r.resamp.bspline, r.in.lidar and r.surf.idw.

The error I keep getting is:

  span style="color:darkgray">ERROR: Raster map
  output2c10671429724b50a23f93a6b91af9fd not found.

The name of the output raster map is different when I change the input file or tool.

First I thought the problem occurred when trying to create a temp output file, so I put it on perm file in a different location but this didn't help.


texts: Want positive Y to be downward (not upward) in OpenLayers3

texts:
I want to use OpenLayers3 simply to draw some shapes that I can inspect by panning and zooming. There would be no map backdrop since the data is not geographic in nature.

The data holds the locations of these shapes using coordinates (X, Y) which have the following convention:

POSITIVE X is EASTWARD from origin
POSITIVE Y is DOWNWARD from origin

However, OpenLayers' default projection has POSITIVE Y going UPWARD.

To try to solve this issue, I am creating a new custom projection and specifying its axisOrientation as 'esu'. However this isn't doing anything. Here is my code...

http://jsfiddle.net/x9ust0nx/ (note there is no map backdrop -- just a grey canvas with 3 red points)

Any ideas?

(I realize I can iterate through my data and negate the Y values, but I want to learn how to do it the right way).


texts: Load vector layer with QgsVectorLayer using MSSQL by select statment

texts:
I want to read a layer from a database. I would like to access the table via a select statment.
In PostgeSQL works both:

vlp1='.....table="mytable" (geom)'
QgsVectorLayer(vlp1, LayerName , "postgres")
vlp2='.....table="(select * from mytable)" (geom)'
QgsVectorLayer(vlp2, LayerName , "postgres")

With MSSQL, only the direct specification of a table works:

vlp1='.....table="mytable" (geom)'
QgsVectorLayer(vlp1, LayerName , "mssql")

.

The select statement does not work:

vlp2='.....table="(select * from mytable)" (geom)'
QgsVectorLayer(vlp2, LayerName , "mssql")

Is the MSSQL Select statement not supported?


texts: Linking polygon features in one layer to polygon features in another, using largest overlap in QGIS?

texts:
I have two polygon layers: the first one contains small land lots (green in the image), the second one contains large zones (the purple line in the image attached separates 2 such zones).

I want to link (the id's of) the large zones in layer 2 to the small parcels in layer one, using the largest overlap as criteria or measure: the small parcel gets the id of the large zone with which it has the largest overlap.

How can this be done in QGIS?


texts: How to update or uninstall new "All-In-One" QGIS macOS packages

texts: Recently new "All-In-One" QGIS packages for macOS have been released. This new packages make installation very simple but inside the .dmg file I have not found any readme file for uninstalling or updating: How can it be done correctly?


texts: How can I render my own tiles using data I get with cartodb's core API library?

texts:
I am trying to use the cartodb.core.js library to render my own tiles. According to what I see, the Tiles library is only available in the cartodb.core.js. And if we use that we don't get the createLayer function, which is fine if I can just get the Tiles to load, and I can render them myself:

var layerData = {
  user_name: 'swc28',
  sublayers: [{
    sql: "SELECT * FROM zoning",
    cartocss: '#layer { marker-fill: #F0F0F0; }'
  }]
};
cartodb.Tiles.getTiles(layerData, function(tiles, err) {
  if(tiles == null) {
    console.log("error: ", err.errors.join('\n'));
  }
  console.log("url template is ", tiles.tiles[0]);
})

But I've never rendered layers like this myself, and I looked at the library without much  success to figure out how to do it. So my question is: How does the CartoDB UI library render the tiles from a template url that is returned form the call above?

The template url looks like this:

http://{s}.ashbu.cartocdn.com/swc28/api/v1/map/fad6aaab9c3d921f76e0e07f3da6071f:1414011659683.92/{z}/{x}/{y}.png?cache_policy=persist


texts: R {spdep} Multiple y observations per region (pooled cross-sectional data) - How to construct weights matrix?

texts:
I'm in a situation where I have a shapefile of N districts (provinces and sub province level) across the globe. I also have a dataset of S survey responses where each individual is tied to a district. I want to account for spatial autocorrelation in a regression but multiple individuals are tied to the same district. For instance, I might have Colorado in the US with 45 respondents.
I want to construct a weights matrix using spdep. I can do this easily with the shapefile data however that matrix is going to only be NxN. The regression functions expect a DxD matrix though (or D long list). I could hard code such a matrix/list from the NxN matrix however I'm wondering if there's a way I could do this in spdep instead, where each individual in Colorado gets the same row from the NxN matrix automatically.
I assume that assigning a polygon (sf or sp packages) to each Survey respondent is not the smart thing to do?


texts: How to convert ArcGIS points into Lat and Long

texts:
I have no experience with GIS data but I have been tasked with taking an export from an ArcGIS database and converting the values into Lat and Long. The client we are helping also has very limited experience with the ArcGIS data they have so it is the blind asking the blind to help them.

Here is a sample of what the client has provided from their ArcGIS database.

1074409.65305250,14531615.53041425,North Garage
1074416.77442934,14531614.91000867,North Garage

What are the values they provided and is it possible to convert these values into Lat/Long?


texts:
I have a linestring in shapely which I'm converting into buffer.

>>> from shapely.geometry import LineString
>>> line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
>>> dilated = line.buffer(0.5)

Now I'm trying to convert the dilated linestring to the original linestring. So I tried doing this:

>>> original = dilated.buffer(0.0)

But it's returning a polygon instead of a linestring. How can I get the original linestring from the polygon in shapely ?